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ABSTRACT 



This thesis presents two methods of ship classification 
with IR images. The first method is the Fourier Coefficient 
method which transform the sample points of a superstructure 
profile of a ship to the spatial frequency components. The 
shape of the coefficient curve can be used to classify the 
type of a ship from 8 different categories. But, the 
differences is so minor that it is difficult to implement a 
computer program to recognize it. The second method is a 
B-spline Coefficient method which uses the uneven spaced 
spline coefficients to find the beginning, the peak, and the 
area of the lumps of a ship for classification. This method 
is better than the Fourier Coefficient method. The study of 
These methods is presented here. 
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I. INTRODUCTION 



It is quite important to classify the ships according to 
their types. Classification can be done in a number of 
different ways. The simplest is the visual method that is 
prone to error. Other methods for classification are the 
Fourier Coefficient method and the B-spline Coefficient 
method. In these methods the necessary information for clas- 
sification can be obtained from the superstructure profile. 

The Fourier Coefficient method samples the superstruc- 
ture profile at every chosen points. The function values at 
the sampling points are transformed into the spatial compo- 
nents. The logarithm of the magnitude of these components 
is plotted and compared with the standard plot to recognize 
the type of the ship. 

In the B-spline Coefficient method, the spline coeffi- 
cients along the X axis and the Y axis are used to recon- 
struct the superstructure profile. The shape of the curve of 
the spline coefficients is in some way similar to the shape 
(the position of the lumps) of the superstructure profile. 
The ship classification may be achieved by recognizing the 
beginning, the peak, and the area of the lumps. 
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II. PREPROCESSING 



Preprocessing is the procedure to obtain the superstruc- 
ture profile of a ship from the IR image. Then, Fourier 
coefficient or spline coefficient methods can be used. The 
details of the preprocessing procedures are a follows. 

A. DATA COLLECTION 

The data consists of the IR image of eight different 
types of ships. 

1. DD - Destroyer; "HALL" class. 

2. Container; The class is unidentified. 

3. Freighter; The class is unidentified. 

4. AOR - Replenishment 'oiler; "WICHITA" class. 

5. LST - Tank landing ships; "NEWPORT" class. 

6. FF - Frigate; "GARCIA" class. 

7. CGN - Guided missile cruiser (Nuclear propulsion); 
"BAINBRIDGE" class. 

8. DDG - Guided missile destroyer; "CHARLES F. ADAMS" 
class . 

These images are taken from an aircraft which is flown 
at a height of 500 feet with a speed of 400 knots toward the 
side of the ship. The aspect angles for these images are 90 
degrees which may be slightly off in some images. The inac- 
curacy of the aspect angle arises from the fact that the 
photos are taken while the aeroplane keeps on moving. All 
data of the images are stored in a digital magnetic tape 
with 256 by 64 bytes per image. Thus, the number of bytes 
required for each image is 16384 and each byte represents 
the intensity of a pixel. For each record of the image, a 
label is coded in the last 8 bytes as follows: 

1. Byte (16377) = Run number in each flight which passes 
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the ship. 

2. Byte (16378) = Video tape time code when the data is 
taken; in minutes. 

3. Byte (16381) = Video tape time code; in seconds. 

4. Byte (16382) = Video tape time code; in thirtieth of 
a second. 

5. Byte (16379) = Range in kilo-feet which is the 
distance measured from the radar it may have an 
error 1 to 2 kilo-feet. 

6. Byte (16380) = Aspect angle; degrees from the bow 
of the ship. 

7. Byte (16383) = Ship class. 

8. Byte (16384) = ID, 1 = for training, 2, 3, 4, 5 = for 

testing . 

The run number and the time code together uniquely 
define a specific image that represents a single TV frame 
with no averaging. In addition, the time interval between 
the end of one image to the beginning of the other is 
approximately 1.5 seconds. Also, there are inherent random 
noise in the record which arises from the photo instrument 
and the process of storing them on to the digital magnetic 
tape . 

B . SOBEL OPERATOR 

The Sobel Operator technique is used to find the edge. 
To determine the edge, the Sobel Operator uses the differ- 
ence of gray levels of the pixels in a 3 by 3 matrix as 
shown in Figure 2.1. 

a , b, c , d, e , f , g , h, and i are the values of the gray 
levels at the position of (x-l,y), (x,y-l), (x-l,y), (x,y), 

(x+l,y), ( x-1 , y+1 ) , ( x, y+1 ) , and (x+l,y+l) respectively. The 

Laplacian estimate is defined as 

2 

d_L ~ f(x,y)-f(x + |, y)-[f(x+ l,y)-f(x+2,y)l 

6x 2 L 1 (2.1) 
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The basis vector for all directions are (a-2b+c), 
(g-2h+i), (a-2d+g), and (c-2f+i). Furthermore, each of the 

basis vector is convolved with the image as follows: 
along the x axis 

d, = [f(x- l,y- 1)+ 2f(x,y-|)+f(x+l,y-l)] . 



-[f(x-l,y+l)+2f(x,y+l)+f(x+l ; y+l)] 



along the y axis 

= [f(x+i,y-i)+ 2 f(x+ u)+f(x+i i y+i)] 
“[Kx-l,y-l)+2f(x-| ; y)+f( x -|,y+|)] 



(2.3) 



Since the magnitude of the resulting vector is the abso- 
lute value of the convolved results, the edge magnitude 
S ( x , y ) [Ref. 1], 



S( x ^y) - (df+djf)/ (2.4) 

Note that the Sobel operator does not use the gray level 
at the position (x,y). The advantage of using a Sobel oper- 
ator over others is that the resulting edge is smoother due 
to a 3 by 3 matrix approach. If we compare the Sobel oper- 
ator with the Laplacian operator, it is seen that the Sobel 
operator using the four basis vector as shown above will 
provide more accurate reading because of noise reduction in 
the original image. Hence, The Sobel operator is often used 
in the preprocessing operation. 
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Figure 2.1 Sobel Operator. 

C. THE USE OF A SOBEL OPERATOR 

When a Sobel operator is used at the edge of the image 
frame, the pixel level which lies out of the frame will be 
set equal to that of the adjacent pixel within the frame. 
The original image of a guided missile cruiser is shown in 
Figure 2.2. Since the result of the Sobel operator is numer- 
ically greater than 8 bits range, we have to rescale the 
result back to 8 bits range. This is achieved by deter- 
mining the maximum and the minimum of the gray level. They 
are then used to rescale the gray level in the results of 
the Sobel operator as shown in Figure 2.3. In some 
instances, the original image is very poor as shown in 
Figure 2.4. Attempts to determine the edge of this image 
failed as shown in Figure 2.5. 



Figure 2.2 A CGN at a Range of 45000 feet. 
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Figure 2.3 Image from a Sobel Operator. 




Figure 2.4 Blured Image of a LST at a Range 67000 feet. 




Figure 2 . 5 



Noisy Image 



from a 



Sobel Operator. 
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D. EDGE THRESHOLD STRATEGIES 

Edge threshold strategies are used to extract the edge 
profiles from the Sobel results. In this case, we use only 
the pronounced value of an edge element at x if g(x) is 
greater than certain threshold value [Ref. 2]. 



' = g(x,y) if q(x y)^. thresho I d 



G(*,y) < 



= 0 



othe rwi se 



(2.5) 



To increase the contrast of the image to a silhouette 
form, G(x,y) is defined as 



G(*,y) 



= 255 



= 0 



if S(x,y) ^threshold 



otherw ije 



( 2 . 6 ) 



The choice of the threshold value is based upon the 
histogram of the edge image as shown in Figure 2.6. The 
estimated critical gray level is chosen so that a majority 
number of the pixels with value between 0 to 255 will fall 
below the critical value. Alternatively, histogram equaliza- 
tion may be used to determine the desired threshold level as 
shown in Figure 2.7. In this case, trial and error method 
was used in conjunction with the above method to obtain the 
threshold so that the correct profile is ascertained. 
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Figure 2.6 Histogram of Figure 2.3. 
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Figure 2.7 Cumulative Distribution of the Histogram 

in Figure 2.6. 
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E. 



HOW TO EXTRACT THE PROFILE 



The edge image of the guided missile cruiser shows 
little variations in the gray level as shown in Figure 2.3. 
These variations are caused by the noise in the original 
image. In this case, the choice of the threshold value is 
based upon both the histogram and the cumulative distribu- 
tion of the edge image so that it contains 90 % of the 
pixels. Therefor, the chosen gray level is 110 and the 
result is shown in Figure 2.8. 




Figure 2.8 Silhouette of a CGN in Figure 2.3. 



F. HOW TO OBTAIN THE SUPERSTRUCTURE 

The original image . of the ship is taken from the aero- 
plane with different displacement from waterline to the 
superstructure in a rough sea, so that the profile of the 
ship with respect to the sea surface varies. Thus, we have 
to eliminate some information in the image of the ship by 
considering the superstructure only. In considering the 
overall ship structure, it is obvious that the largest 
distance is between the bow and stern span. Therefore, the 
bow and stern points are located first in the program as 
shown in Appendix A. Then we consider the slope of the bow 
and stern of the ship and set all the gray values below the 
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slope line to zero. Hence, it results in a superstructure of 
the ship as shown in Figure 2.9. However, in some images, 
there is a lot of noise in the background which cause diffi- 
culty in locating the bow and stern. Under these circum- 
stances, the dimensions of a ship are estimated by trial and 
error method. Then the gray level which lies outside is set 
to zero and an estimate of the bow and stern slope can be 
made . 




Figure 2.9 Superstructure Profile of Figure 2.8. 



G. CONTOUR FOLLOWING 

The noise in the background of the original image, 
yields wide edge structure. In considering this factor, the 
contour following procedure tracking the inner part of the 
image in Figure 2.9 gives the superstructure profile. The 
objective of this contour following is to describe the bow 
and stern points of the ship, the direction, and the posi- 
tion of the edge of the superstructure. The contour tracing 
is done in the counter-clockwise direction which compares 
pixel value of 0 or 255 in a 3 by 3 matrix in the following 
manner . 
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Starting from the leftmost point in the superstructure 
image in Figure 2.9 the contour profile tracing is accom- 
plished by examining the neighbors of a 3 by 3 kernel 
located at the curser position as shown in Figure 2.10. The 
curser is moved along the profile. All successive positions 
of the curser constitute the contour profile of the ship. 
The testing procedure is explained below. 

1. Initialize the curser position to the beginning of 
the thresholded image with the gray level of 255 and 
the estimate direction of the next position. 

2. Check for reaching the end of the profile, if it is 
at the ending of the profile then stop, if not go 
to 3 . 

3. Check for the estimate to see whether it is in the 
direction of. North, East, South, or West. If the 
direction it is North then go to 4, if it is East 
then go to 5, if it is South then go to 6, and if it 
is West then go to 7. 

4. Determine the next position with the gray value of 
255 in the counter-clockwise direction as shown in 
Figure 2.11. Store the position found and move the 
curser to that position. Update the estimate 
direction to the one last found; then go to 2. 

5. The procedure is the same as that in step 4 except 
search pattern is shown in Figure 2.12. 

6. The procedure is the same as that in step 4 except 
search pattern is shown in Figure 2.13. 

7. The procedure is the same as that in step 4 except 
search pattern is shown in Figure 2.14. 

The flow chart of the procedure is shown in Figure 
2.15, and the detail of each procedure are included in 
Appendix B. The testing procedure have to be performed in 
such a maner that the resulting contour is a good represen- 
tation of the superstructure line. The result of the contour 
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image is shown in Figure 2.16. If the superstructure in 
Figure 2.9 has wide edge, we can not find the contour image. 
We have to use the additional step which is called the 
"Closing" operation. 




Figure 2.10 The 3 by 3 Kernel. 
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Figure 2.11 The Step Checks in the North Direction. 
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Figure 2.12 The Step Checks in the East Direction. 
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Figure 2.13 The Step Checks in the South Direction. 
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Figure 2.14 The Step Checks in the West Direction. 
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initialize the beginning position 
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Figure 2.15 Flow Chart of Contour Following. 




Figure 2.16 Contour Image of a CGN at 
a Range of 45000 feet. 



H. CLOSING OPERATION 

Some results from the superstructure extraction process 
are discontinuous because the gray level of those areas of 
the structure is less than the threshold value. If we 
decrease the threshold value, the details of the superstruc- 
ture are effected. It is necessary to use the "Closing" 
operation which consists of the "dilation" process to smooth 
the superstructure profile used the "Erosion" process. All 
direction are dilated similarly which causes an smoothing 
effect on the edges. The superstructure increases in total 
area. Then, use the "Erosion" process to shrink (subtract) 
the dilated part in all directions, thus obtain the smoothed 
superstructure with the appropriate size. The Closing 
process employs the dilation process which yields an output 
image from an input image. The Dilation results is shown in 
Figure 2.17. 
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Figure 2.17 The Process of Dilation. 

We examine the gray value of each pixel in the original 
image. Begining from the pixel at the first row and the 
first column. The procedures are the following. 

1. Considering one pixel in the original image with the 
kernel (B) centered there. If at least one pixel in 
the kernel has a value of 255, we let the gray level 
of the output image at the center of the kernel be 
255. 

2. We shift the center of the kernel 1 column to the 
right. Then following the same procedure as in step 1 
until the last column is reached. 

3. We shift the position of the kernel to the next row 
and starting from the first column. Then, following 
the same procedure as in step 1 and step 2 until the 
last row and the last column is reached. 

The result obtained from the dilation process is an 
image with enlarged structure. The second procedure in the 
Closing operation is the Erosion process. The Erosion 
process perform the same procedures as the dilation process 
except for step 1. If every pixel in the kernel are 255, we 
let the gray level in the new image at the center of the 
kernel be 255. Otherwise, it will be 0. The output image 
obtained will have smooth edge with minor change occurring 
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in the edge detail as shown in Figure 2.18. In this case, we 
use an kernel of size 3 by 3. I f we increase the size of the 
kernel, the details of the image are decreased. 




Figure 2.18 The Profile after Dilation and Erosion 

(Closing Process). 



PROFILE ROTATION 



Often in the contour image, the first and the last point 
of the superstructure are not at the same horizontal level. 
We have to rotate the contour image by setting the two 
points to the same horizontal level. How to rotate it from 
the Y, X axis to the Y', X' axis is shown in Figure 2.19. 

If the angle value Q is positive, it will be 



"x 




cosQ 


-sin^l 




‘X 


Y'_ 




sinQ 


cosQ _ 




_Y_ 



If the angle value Q is negative, it will be 
cos Q sin^" 
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J 




__ 1 

X 

1 





cos 



e. 





*x' 




„Y. 



For some of the contour profile, part of the profile 
after rotation will be out of the image frame. Then we have 
to shift this contour down by 20 pixels position. The 
rotated profile is shown in Figure 2.20. 



29 




'I 

\ 


/ 

t 

\ 

^T\ ^ 

1 \ 

, v' 







Figure 2.19 Rotation Process. 




Figure 



2.20 Rotated Profile of 
a Range of 45000 feet. 



a CGN at 
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III. FOURIER COEFFICIENT METHOD 



We have obtained the ship profile in the previous 
chapter. To extract features out of this profile for clas- 
sification purposes, we will use the Fourier Transform 
method. 

The Fourier transform of the ship profile showed that 
the transform coefficients depend upon the ship's dimen- 
sions, its superstructure, and the distance between the 
camera and the ship. If the profile f(x) is a discrete 
function with 128 sample points. The discrete Fourier 
transform can be written as 

F(u) =1£ f(x)e xp(-j27Tu x/N) (31) 

For u, x = 0, 1, 2, , N-l. N is the total number of 

samples . 

If the direct calculation of the discrete Fourier trans- 
form is chosen, the number of complex multiplication and 
addition will be equal to N ; i.e. to obtain F(0) would 
require complex multiplication and addition N times. In 
order to reduce the computation, we use the fast Fourier 
transform algorithm. Thus, equation 3.1 [Ref. 2] can be 
separated into F (Ven (u) and F odd (u) 

F„,»=ip(2x)^ X 

x=0 

F 0 „(u)=iyf(2x + I)w: x (3.3) 

Ml—, 

x=o 
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W = e x p(- j 27T/M) 

M =_N 
2 

F(u)=i[F fKf „(u) + F^ tf(u )W 2 U M ] 
F(u + M)=l[ F „ f »-F otf ,(u)^J 

Using this method, we have reduced the total number of 
complex multiplication and addition to Nlog N. In this case, 
we have N = 128, thus, the total number of complex multipli- 
cation and addition will be 896. 

First, we divide the rotated profile image into 128 
divisions. Since the distance of each division is equal to 
the amount of the pixel between the bow and the stern of the 
ship divided by 128, we use the distance perpendicular from 
the horizontal line between bow and stern to the highest 
point of the superstructure as the sampled values. The 
result of the fast Fourier transform is a complex number. 
Then we use 



(3.4) 

(3.5) 

(3.6) 

(3.7) 



M(u) - lo g[l + G(u)] (3.8) 

G(u) is the magnitude of F(u). A value of 1 is add to the 
magnitude to avoid negative logarithm result, the results 
obtained are as shown in Table I and 



1 . 


Figure 


3.1 - Destroyer 


2. 


Figure 


3.2 - Container 


3 . 


Figure 


3.3 - Freighter 


4. 


Figure 


3.4 - Replenishment oiler 


5. 


Figure 


3.5 - Tank landing ship 
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6. Figure 3.6 - Frigate 

7. Figure 3.7 - Guided missile cruiser 

8. Figure 3.8 - Guided missile destroyer 

On the results minor difference of the shape of the 
Fourier coefficients can be noticed visually. But it is 
difficult to implement a program to detect these minor 
difference in shape. Therefore, a Second method was used to 
handle this problem. 
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Figure 3.1 Logarithmic Magnitude of the Fourier Transform of a DD. 
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Figure 3.2 Logarithmic Magnitude of the Fourier Transform of a Container. 
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Figure 3.3 Logarithmic Magnitude of the Fourier Transform of a Freighter. 
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Figure 3.4 Logarithmic Magnitude of the Fourier Transform of a AOR. 




6 



38 



Figure 3.5 Logarithmic Magnitude of the Fourier Transform of a LST. 
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Figure 3.6 Logarithmic Magnitude of the Fourier Transform of a FF. 
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Figure 3.7 Logarithmic Magnitude of the Fourier Transform of a CGN. 
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Figure 3.8 Logarithmic Magnitude of the Fourier Transform of a DDG. 



TABLE I 

Logarithmic Magnitude of Ships 



F(N) 

(IE-1) 


DD 

(IE-1) 


Cont. | 
(IE-1) | 


Frig . 
(IE-1) 


| AOR | 
1 (IE-1) | 


L3X 

(IE-1) 


1 ff 

1 ( IE- 1 ) 


CCN 

(IE-1) 


DDG 

(IE-1) 


DC 


1 . El 


1.E1 | 


1 . El 


1 1 . El | 


1 . El 


| 1.E1 


1 .El 


1 .El 


F(l) 


6.74 


7.23 | 


3.66 


1 4.28 | 


6.81 


| 6.97 


7.78 


8.10 


F(2) 


2.17 


7.18 | 


5.17 


| 8.04 | 


2.64 


| 5.67 


6.37 


5.04 


F(3) 


3.21 


6.45 | 


7.88 


| 3.84 | 


3.35 


| 5.64 


2.38 


5.40 


F( 4 ) 


3 . 17 


6.18 | 


3.76 


1 5.28 | 


1.28 


| 4.43 


7.6E-1 


3.26 


F(5) 


2.20 


5.09 | 


5.31 


1 2.45 | 


2.49 


| 2.40 


3.53 


4.81 


F(6) 


1.88 


4.57 | 


3.17 


1 2.66 | 


1.36 


| 1.52 


3 . 48 


6.02 


F ( 7 ) 


1.10 


2.03 | 


3.90 


| 2.28 | 


3 . 6E- 1 


| 8 . IE- 1 


3 . 44 


4.50 


F(8) 


1.50 


2.65 | 


4.58 


1 1-07 | 


1.74 


| 1.76 


3.35 


3 . 14 


F(9) 


2.07 


1.83 | 


1.67 


1 1-32 | 


2.14 


I9.9E-1 


3 . 40 


2 . 46 


F( 10 ) 


1.13 


1.57 | 


2.07 


| 3 . 39E1 | 


1.19 


| 1.13 


3.24 


1.05 


F ( 1 1 ) 


4.6E-1 


6. 0S-1I 


4.53 


1 1-03 | 


2.46 


| 1.14 


2.62 


1.23 


F( 12 ) 


1 . 6E-1 


5.6E-1 | 


2.23 


1 1-07 | 


1.73 


I9.8E-1 


1.80 


2.09 


F( 13 ) 


5.4E-1 


1.62 l 


2.09 


| 5.0E-1 | 


1.78 


| 3 . 5E-1 


a. 5e-i 


1.23 


F( 1 4 ) 


9 . 6E-1 


1.48 | 


2.09 


I9.8E-1 | 


1.87 


| 7 . 4E- 1 


6.7E-1 


1.62 


F( 15 ) 


4. IE- 1 


1.25 | 


1.07 


| 7.CE-1 | 


1.14 


| 4. IE- 1 


1.19 


1.38 


F(16) 


2.9E-1 


1.77 | 


1.30 


1 1.05 | 


7 . IE- 1 


| 1.36 


1.25 


1.30 



DD = Destroyer at range 77000 feet. 

Cont = Cointainer at range 28000 feet. 

Frig= Freighter at range 40000 feet. 

A0R = Replenishment oiler at range 78000 feet. 

LST = Tank landing ship at range 51000 feet. 

FF = Frigate at range 49000 feet. 

CGN = Guided missile cruiser at range 45000 feet. 
DDG = Guided missile drestroyer at range 41000 feet 
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A. VARIATION OF SHIP SUPERSTRUCTURE WITH RANGE 



One of the practical problem of using the ship profile 
is that it is sensitive to range variations. Close ship 
profile has more details than the far away ship profile. 
The dependency of the geometric size of the profile with the 
range is discussed in this section. 

Assume that the object is centered on the camera axis. 
The field of view of the camera, the number of the pixel in 
the image, the size of the image, and the size of the object 
are known. Our problem is to determine the distance between 
the camera and the object. 

1 . Background 

The camera system is similar to a human eyes system. 
The light reflected from the object goes into the eyes. The 
image of the object falls on the retina, the signal is sent 
to the brain in some electrical from, and the brain changes 
it to a from that human can perceive. But, in the camera 
system film or senser are used to pick up the image. The 
function of a camera is shown in Figure 3.9 




Figure 3.9 The System of the Camera. 



43 




Figure 3.10 Range Calculation. 

The distance between the lens and the image plane 
can be adjusted in order to have a clear image on the film. 
The camera has a field of view (FOV) angle as shown in 
Figure 3.9. The object has to be in the field of view of 
the camera. The determination of the distance is shown in 
Figure 3 . 10 

For simplicity of calculation the inside of the 
camera is flipped to the same side as the object as shown in 
Figure 3.10. when the angle of the FOV isO(, 1/2 is the 
half of the full image size, D/2 is the half of the dimen- 
sion of the object, X is the distance from the lens to the 
image, and R is the distance from the lens to the object. 

Assume that the distance 1/2, distance Im/2, and 
angleQ^/2 are known. Then, the distance R can be determined 
by 



X - JjiLta na 
2 2 



ianQ' , = IJ_ 
2 2 X 



(3.10) 
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(3.11) 



R = D I 
2 tana 
2 

Assume that the angle resolution of the pixel of the 
field of view of the camera is equal to 0.2E-3 radian per 
pixel. The number of pixel of the frame is equal to 256. 
The size of an image is I pixels. The dimension of the 
object D in feet is known. The field of view in angle is 



a =(0.2E- 3)256 



X =.256. c[ 

2 tanQL (3.13) 

2 



where d is in unit of pixel. 

tan a; = I d 

2 71 ( 3 - 14 ) 



tana - I tana 

2 256 2 (3.15) 



R =JL L 

2 tong 1 
2 



(3.16) 



R = 1 28D 
Ita n 0.0256 



(3.17) 



When the length D of the object is known, from equa- 
tion 3.17 we can estimate the distance from lens to the 
object and is shown in Table II 
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TABLE II 

Range Estimation 



CLASS 


I ( pixel ) 
MEASURE 


1 D (ft) 
| KNOWN 


1 R (kft) | 
1 1 
| CALCULATE | 
1 1 


i 

R(kft) | (R-R)IOO 

1 

RADAR DIS | R 

1 


DDl 


96 


| 418 


1 1 

I 21.766 | 

1 1 


77 


i 

I 71.71 
1 


DD2 


80 


418 


1 1 

| 26.119 | 

i i 


85 


1 

| 69.27 

| 


A0R1 


107 


| 659 


I 1 

| 30.787 | 

1 


78 


1 

l 60.53 

i 


A0R2 


96 


| 659 


| 34.315 | 

i i 


85 


1 

| 59.63 

i 


LST1 


176 


| 522.3 


1 1 

| 14.834 | 

i i 


51 


1 

| 70.91 

i 


LST2 


134 


| 522.3 


1 1 

| 19.484 | 

i i 


57 


| 65.82 

1 


CCN1 


147 


| 565 


1 1 

| 19.213 | 

i i 


45 


| 57.31 

i 


CGN2 


119 


| 565 


1 1 

| 23.734 | 

i i 


55 


1 

| 56.85 

I 


DDG1 


126 


| 437 


1 1 

| 17.337 | 

i i 


47 


| 63.11 

i 


DDG2 


90 


| 437 


1 1 
| 24.247 | 


64 


1 

| 62.08 



DDl , DD2 = Destroyer at range 79000 and 83000 feet. 

A0R1,A0R2 = Replenishment oiler at range 78000 and 88000 feet. 
LST1 / LST2 = Tank landing ship at range 51000 and 62000 feet. 
CCN1,CGN2 = Cuided missile cruiser at range 45000 and 64000 feet. 
DDG1,DDG2 = Guided missile destroyer at range 41000 and 64000 feet. 



The error distance between the estimated distance 
and calculated distance in percentage is ( ( R - R)/R) 100. 

2. Remark 



Calculated distance error in R may come either from 
the pixel measurement in an image or from the angular reso- 
lution estimation. The distance that is stored in the image 
label has an error of about 1 to 2 kilo-feet. If the angle 
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of the field of view is accurate, we can estimate the number 
of pixel in an object image. Then we can classify the object 
by comparing the number of pixel of the original image with 
that of the test image. Some known system parameters can 
help to determine the range of the target ship. The problem 
is that existing errors in the system parameter causes in 
the larger errors in the estimated range. Consequently, they 
are not very useful in classifying the ships. 
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IV. B-SPLINE COEFFICIENT METHOD 



Using B-spline coefficient is another method to describe 
a ship profile. This method uses the B-spline coefficients 
to determine the beginning, peak, and area of lumps which 
contain significant information about the type of the ships. 
The comparisons of the knot position (in parametric value) 
from the midships to the peak or beginning of the lump, can 
be very helpful. Different ships will have different lump 
positions and sizes. 

A . BACKGROUND 

A spline function is a piecewise polynomial used to 
interpolate points. This kind of curve is smooth and the 
discontinuities in its k-th derivative is as small as 
possible. In this case. Cubic spline was used, where the 
first and second derivatives for any set of interpolating 
points are continuous, while the third derivative may be 
discontinuous. The reason for using Cubic spline is to keep 
the third derivative discontinuity as small as possible and 
the curve as smooth as possible. The B-spline calculation 
procedure is also very stable. In our case the order of 
spline used is 3, while the 1-st and 2-nd derivative are 
continuous. The choice of 4-th order spline function is due 
to the fact that it is generally sufficient for most ship 
profile curves. Discontinuity, in a sense, may be stated as 
the jumps of the third derivative, which is the means to 
control smoothness of the connecting pieces. 
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B. B-SPLINE APPROXIMATION WITH FREE KNOT 



As mention earlier, the B-spline function is used in the 
spline approximation with free knot. The free knots is some- 
times called uneven or irregular knots; that is no fixed 
number of knots are used. Furthermore, the spline position 
need not be on the original curve so as to minimize the 
number of knot position while preserving most of the detail's 
of the original curve. 

First, minimize the value of the lack of the smoothness 
Tj(c) defined as [Ref. 3] 

77(c) = Z ( X o. C . ) 

J=l <=-*•' (4.1) 



where Ci is a coefficient of spline at the knot position, 
a ;j is defined as follows 

= (4 . 2) 

where M ; - ^ is the normalize B-spline function and defined 
as 

«;,*+/(*) = (W *+r t;)A* + '(t, , • • • ,t, + 1+ , ) 9,<t; x) 

(4.3) 

and G* ( t; x) , t are defined as follows 



(=((-*)*= (t-x)* 




if x 



if t< x 



(4.4) 
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& L f (Z;,Z ,>/ , . . . ,Z )f(t) stands for the 1-th divided 
difference of the function f (t) on the point Zj , . . . , Z ( - + ^ 
where t is the position values of knots in term of Z param- 
eter. Z parameter is defined as follows 

Z(L) =Z(I-I)+ [(X(l)-X(l-l)f+(Y(l)-Y(I-l)) 2 ]^ (4.5) 



Z(0) = 0 



(4.6) 



where I is the number of the sampling point and 
1=1 , 2 , 3 , . . . , m. 

Second, the smoothing is subjected to a constraint 



5(C) < S 



(4.7) 



where S is the smoothing factor, 0\C) is the weighted sum of 
the square residuals defined as 

(4 8, 

j-l '=-k 

Xy , Yy are the values of X and Y at Z parameter of the 
sampling points; W- is a weighting factor for all sampling 
points [Ref. 3] defined as 



W‘ 




(4.9) 



where W- is an estimate of the standard deviation of Y- . 

/ I' ■— J 

Then, the value of S is in the range m+V2m. If nothing is 
known about the statistical standard deviation of Yy , each 
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W can be equal to one and the S can be determined by the 
trial and error method. 

1. Interpolation and Approximation Using B- spline 

Function 

There is a difference between the interpolation and 
the approximation. In interpolation, the number of knot 

positions required are equal to the number of sampling 
points and the value of S in eq(4.7) is small. In this case, 
the computation time required will increase tremendously. 
Whereas, in approximation, it is not of great concern that 
the approximated value at a position be the same value of 
the sampling point, and most of the information still remain 
in the original curve. Thus, decreasing the value of S will 
result in the large number of the knot positions and the 
final curve obtained will be similar to the original curve. 

The justification for using approximation rather 
than interpolation approach is that though the resulting 
curve may not be the best fitting curve, it is smooth and 
close enough to the original. The approximation spline func- 
tion has less number of knots than the number of samples, 
which reduces the total processing time. In approximation 
approach, when the appropriate choice of S value is made, it 
will, in turn, generate sufficient number of knot positions 
required to provide a close approximation to the original 
curve. The ratio of sampling point to spline coefficient is 
10 to 1 as shown in Figure 4.1. In Figure 4.1, a guided 

missile cruiser ship at a distance of 45000 feet, with 290 
sampling points and the associated B-spline knots are shown. 
The original samples are plotted in solid curve. The knot 
position are show by small circles in Figure 4.1. After 
B-spline approximation, the number of the knots are reduced 
to 33. Hence, this technique is essentially a kind of data 
compression technique. 
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Figure 4.1 
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Plot of the Original Data and the Approximate 
with Free Knot. 

C. TO DETERMINE THE KNOT POSITION AND THE B- SPLINE 
COEFFICIENTS 

The approximate smoothing factor in eq(4.7) is to be 
calculated before using the subroutine "PARAM" [Appendix C]. 
In the subroutine, there is a check on S factor every time 
it is run. If the S input values do not satisfy or meet the 
criteria, the program will return with a error code IER. 

There is a parameter NEST which is set to a constant 
value. This value determine the dimension of the knot posi- 
tions which relates to the array T(NEST), Cx(l..NEST) and 
Cy(l..NEST). The NEST is an overestimate of the dimension of 
the arrays set by the user. The limitation on the value of 
NEST for subroutine "PARAM" are as follows [Appendix C] 

1. 2k+2 NEST M+.k+l 

2. Typically, value of NEST M/2 

where M is the number of the total sampling points and k=3 
is the highest order in the B-spline function. 

The subroutine "PARAM" produces several outputs as, N 
the total number of knot positions, T(N) the array of the 
value of knot positions in Z parameter, Cx(N) and Cy(N), 
B-spline coefficient of X functions and Y functions for knot 
positions defined in array T(N). 
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1 . 



Limitation on B- spline Coefficient Determination 



If the value of NEST is too small, the user will 
receive error code IER and the number of knot positions will 
appear to be dense in the first part of the curve. While 
sparse in the other part. The example of an incorrect selec- 
tion -of NEST is shown in Figure 4.2. 




Another problem which causes difficulty in program- 
ming is that the main program is in PASCAL while the subrou- 
tine is in FORTRAN. As for the PARAM program in FORTRAN, the 
array index value starts from 1, while , for the user PASCAL 
program it starts from 0. Therefore, in linking the main 
program to the subroutine, one has to keep in mind of the 
difference. In addition, the FORTRAN programming logic 
structure is so complicated that, when an error occurs, it 
is very difficult to debug and locate the error. 

The values of Cx and Cy depend upon uneven knot 
positions, and they contribute controlling effect to the 
reconstructed curve which would be both smooth and close to 
the original curve. In running the B-spline approximation 
program for the first time, the values of Cx and Cy of the 
last 4 knots at the right end are zeros. However, in running 
it again, increase the value of S did not yield zero values 
of Cx and Cy at those points. This, nevertheless, has no 
effect on the reconstruction of the curve. 
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Plots of the reconstructed profile of three ships of 
the CGN class using different value of S, result in 
different number of knot positions and are shown for compar- 
ison in Figure 4.3, 4.4 and Figure 4.5. The original 
samples are plotted in solid curve, the reconstructed curve 
is plotted in dashed line in Figure 4.3. 
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2. Selection of Smoothing Factor Value ( S ) 

It is found by trials and errors using the computer 
program PARAM that in order to retain maximum information of 
the profile curve, the number of knot positions required 
should be in the range of 25 to 35. The number of knot 
positions depends upon the value of S which must be set 
accordingly. The appropriate choice of S to satisfy the 
condition stated above, is important. For the class of a 
guided missile crui ser ( CGN ) , three ship images at 3 
different ranges were selected. Then, run the appropriate 
program for various S factors to see how the number of 
knot(N) will vary. The results are shown in Figure 4.6. 
Plots of N vs S in Figure 4.7 through Figure 4.13 show that, 
in most cases, the value of N decrease quite rapidly when 
the value of S is in the range of 0 to 100, and gradually 
for S factor in the range of 100 to 200, thereafter, the 
value of N decreases very little. Obviously, the curve seems 
to decay exponentially. Furthermore, for some classes of 
ships changes are more pronounced than the others which is 
probably due to the actual number of knots present in the 
profile. For guided missile cruiser ship(CGN) with 2 lumps, 
the number of knot positions required can be 33 as shown in 
Figure 4.6 We select the factor S to be about 100. 

The selection of the factor S depends upon the 
number of the original sampling points. If the factor S is 
small, large number of knots are needed. When the factor S 
is large, small number of knots are needed. When the number 
of knot and the Cx and Cy coefficients are small, the 
B-spline coefficients Cx and Cy obtained can not be used to 
reconstruct the curve close to the original profile. 
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Figure 4.6 Plot N vs S for a CGN. 
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Figure 4.8 Plot N vs S for a Container. 
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Figure 4.10 Plot N vs S for a AOR. 




O'OS O'Ob O'OC O’OS 0 ‘ 0 1 0*0 

(N) lONM JO JjgwnN 



64 



50.0 100.0 150.0 200.0 250.0 300.0 350.0 

SMOOTHING FRCTOR (S) 



o 




(NJ 10N>H JO JjgwriN 



65 



Figure 4.12 Plot N vs S for a FF. 
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3. The Output Cx,C_y and Knots Profiles 



From the previous study the factor S is selected for 
each of the class of ship as follows 

1. Destroyer (DD) , S = 20.0 

2. Container, S = 100.0 

3. Freighter, S = 100.0 

4. Replenishment oiler, S = 20.0 

5. Tank landing ship(LST), S = 25.0 

6. Frigate ( FF ) , S = 100.0 

7. Guided missile cruiser (CGN) , S = 125.0 

8. Guided missile destroyer (DDG) , S = 125.0 

The plot of the B-spline coef f icients , Cx and Cy, and 
X,Y at the positions of the sampling points vs the Z param- 
eter are shown in Figure 4.14 through Figure 4.21. 
Observation and comparisons of the curves show that the 
values of Cx and Cy exhibit' changes similar to that of X and 
Y except that the variation of values leads that of the X 
and Y. This is due to the fact that Cx and Cy have to act 
as the controlling factor for the reconstructed B-spline 
curve to get the result close to the original curve. 

Examination of the plots of the results of X and Y 
show that when X is increasing monotonically , Y is almost 
constant; but when X is almost constant, Y is increases or 
decreases. This behavior relating to profile reconstruction 
may be explained as follows. For a ship profile when X is 
increasing and Y not increase too much, this may be inter- 
preted as an almost leveled profile. When X is almost 

constant, and Y may be increasing or decreasing, it may be 
interpreted as the beginning or the ending of the lump. 
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Figure 4.14 Plot X,Y,Cx, and Cy vs Z for a FF at a Range of 43 K-ft. 
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Figure 4.15 Plot X,Y,Cx, and Cy vs 2 for a LST at a Range of 51 K-ft. 
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Figure 4.16 Plot X,Y,Cx, and Cy vs Z for a CGN at a Range of 45 K-ft. 
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Figure 4.17 Plot X,Y,Cx, and Cy vs 2 for a DDG at a Range of 41 K-ft. 
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Figure 4.18 Plot X,Y,Cx, and Cy vs Z for a DD at a Range of 77 K-ft. 
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Figure 4.19 Plot X,Y,Cx, and Cy vs 2 for a AOR at a Range of 78 K-ft. 















N 












i 


1 




i 


1 1 


! 






1 


1 II 


i 














_a 


L 




























1 


Li 


k 


i 














\\\\ 


























1 


P 6 


















iJs 


r 














i 


1 








f i 


1 <J 


f 


















1 lo 










! 




i 


1 1 


1 1 


1 ii 


l 1 




















! c 


\ 














III! 1 


1 9 




1 


1 
















1 ! * 


k 














1 1 l 1 


j 


i 


\\ 




















i m. 
















1 


i 


1 « 


\ 1 


i 




| 














1 I 




\ 






1 1 
















II 






1 












1 


1 ' 




\ 


1 1 1 


1 


I 


1 1 




i 1 


p ii 


















i 




k) 




\l 


1 






1 






i i 














| 


i 




I-- 




i i 


K 








i 


IIM 


















i 


i 










O.l 
1 O' 






i 






1 lo 


i 














Mil 




i 




1 !U 




i i 


1 1/* 








4* 










1 1 




i 


1 


! J 


14 1 1 l l 


1 i M 




























1 M 




l 






! f. 




























i 






p i 


1 




1 1 Q l N 


r 


i 












I 




i i 


1 














M 


i 






| ; 




i 


i i 


















i 

! 1 












'—O'. ^ 

1 


kl i 


MM 


i ii 


| 


































A 


K 




I 






r 




































>• , 


D 1 1 


l\ 








i 




























I 












io 17 






ITT 






































| 






V 




> 


i 


£> 




















i 






















o 


) 




k .o 


































j 




i 


1 






o 


\ 


tf?l 










a 

e — 






























i 








f\ 


M 


1 , 








t- 








1 
















i 


I 
















H 


M 








<T 








1 


1 


i 


1 | 


| 




















1 






\ 








tx 

K 




1 


| 


1 


I 1 


i 1 


1 1 












1 










□ 












— Cfc 

c 

— Ll 



NOIlISQd g § o 



J 



74 



Figure 4.20 Plot X,Y,Cx, and Cy vs Z for a Freighter at a Range of 40 K-ft'. 
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D. 



SHIP DESCRIPTION 



The shape of the ships depend upon the shapes and posi- 
tions of the lumps. Characteristics of the lumps for 

different type of ships are as follows: 

1. Frigate - The beginning of the lump is at 1/6 of the 
ship length to the left side of the midships and 
the peak of the lump is at the mast which is located 
at the midships. The termination of the lump 

is approximate 1/3 of the ship length from the 

midships to right side. In addition, the average 
height of the lump is a little higher than the 

level between the bow and the stern, and its size 
is 1/2 of the ship length as shown in Figure 4.22. 




Figure 4.22 Frigate. 

2. Container - The beginning of the lump is at 1/3 of 
the ship length to the right side of the midships 
while the lump is high and terminate at the stern. 
The lump appears to be in a rectangular shape with 
small crane. 

3. Tank landing ship(LST) - The beginning of the lump is 
at 1/4 of the ship length from the midships to the 
left side; the height of the lump is higher than the 
level between the bow and the stern by a small 
margin, while its highest point is located at 
approximate 1/6 of the ship length from the midships 
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to the left side with small difference from the 
average high as shown in Figure 4.23. 



4. Guided missile destroyer - the beginning of the lump 
starts at approximately 1/4 of the ship length from 
the midships to the left side with its highest point 
at the mast. After this the slope will decline in a 
very rapid fashion following a noticable deck 
distance, then the beginning of the second lump 
occurs due to the redome presence. Therefore, the 
lump will be narrow with great height, and will 
terminate at approximatly 1/6 of the ship length from 
the midships to the right side. Furthermore, there 
is a little lump near the stern, this 
distinguish destroyer of the same size as shown in 
Figure 4.24. 

5. Destroyer-lump characteristic will be similar to that 
of guided missile destroyer except for the small lump 
as in Figure 4.25. 

6. Guided missile cruiser - The beginning of the first 
lump is at 1/12 of the ship length from the midships 
to the left side with highest point at the mast. The 
lump size is large both in length and height and its 
height decreases to the point which is a little above 
the level between the bow and the stern. Therefore, 
the second lump begins with almost the same size as 




NEWPORT. ' NEWPORT" Cl«i LST 



USA :20t 



Figure 4.23 Tank Landing Ship(LST). 
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Figure 4.24 Guided Missile Destroyer (DDG) . 




Figure 4.25 Destroyer. 

the first one. The second lump ends at approximately 
1/4 of the ship length from the midships to the right 
side. Furthermore, the distance between the peak of 
both lumps is less than that of the Guided missile 
destroyer shown .in Figure 4.26. 




Figure 4.26 Guided Missile Crui ser ( CGN ) . 



78 



7. Replenishment oiler (AOR) - There are 2 little lumps 
with the first one starting at 4/5 of the ship length 
from the midships to the left side exhibiting rapid 
decrease in height to the level between the bow and 
the stern. The second lump starts at 1/4 of the ship 
length from the midships to the right side with slow 
decline to the point near the stern. 

8. Freighter-with 3 noticeable lumps, the first two 
represent lifting crane with the first high lump 
approximately 1/6 of the ship length from the 
midships to the left side but narrow; meanwhile, the 
second lump is located at approximately the 
midships with the same size as the first one. The 
beginning of the third is at approximately 1/3 of 
the ship length from the midships to the right 
side with large size. It is higher than the first 
two and with gradual decline to the stern. 

From the lump characteristic of different type of ships, 
we could distinguish the type of a ship based on the rotated 
superstructure . 

E. SHIP CLASSIFICATION 

An important aspect in categorizing types of ships is to 
recognize the lumps above the main deck. Therefore, it is 
useful to plot the positions of X and Y vs the Z parameter 
where the Z parameter can be determined from eq(4.5) where 
Z(I) is an array (1..M) as shown in Figure 4.27 where X is 
ploted in solid line and Y is ploted in dash line. Then, 
plot the B-spline coefficients Cx and Cy vs T(N); the knots 
position in Z where N is the total number of knots as shown 
in Figure 4.28. The Cx is ploted in solid line. The Cy is 
ploted in dash line. The Values of Cy will be close to the 
curve of Y while the values of Cx for the same knot posi- 
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tions as Cy is exhibiting monotonic increasing trend. When 
the difference in Cx is decreasing or zero, the value of Cy 
will have a pronounced change where the beginning or the 
ending points of a lump can be detected. The value of knot 
position(T) at those points related to the sizes of the 
Cx,Cy can be determined. Finally, with the values of Cx and 
Cy at those points known, the area of the lumps can also be 
determined . 

Hence, from the ship's characteristics and information 
derived from the above procedure, classification of ships 
can be made by considering 

First, the number of lumps detected, 1, 2, or 3. 

1. The number of lump=l: Frigate, Tank landing ship 

2. The number of lump=2 : Destroyer, Guided missile 
cruiser. Guide missile Destroyer, and Replenish- 
ment oiler (AOR) 

3. The number of lump=3 : Freighter and Container 
Second, the position of a lump relative the midships is 

measured. This quantity is scaled by the total length of 
the ship. This scaled quantity will be invariant with 
respect to the different ship sizes at different ranges 

Third, the area of the lump is normalized to the ship 
length squared. 

1 . Lump Detection 

As shown in the plot of X, Y, Cx, and Cy vs Z param- 
eter in Figure 4.15, when the difference (ACx) between 
successive value of Cy varies from increasing to decreasing, 
and then, to increasing again, Cy exhibits noticable varia- 
tion. From observation of Cx values, it is seen that the 
difference (ACx) always has variation in the same sequence 
as stated above. Therefore, only Cy values are taken into 
consideration in the program that detects lump. 
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There are two distinct procedures used in dealing 
with two different kinds of lumps big or small. Thus, it is 
necessary to distinguish in the first place whether the lump 
is big or small. For the big lump, the differences of Cy is 
positive for all initial four knots. It continues to the 
peak and decreases toward the ending of the lump. For the 
small lump, the difference of Cy is positive for the first 
two knots and stay constant or positive for the third knot, 
but the difference of Cy for the forth knot is negative, 
thereafter, the program proceeds to establish the following 
values for each lump detected: 

First, the knot positions (Z value) at 

1. the beginning of the lump 

2. the ending of the lump 
Second, the Knot number for 

1. the beginning of the lump 

2. the ending of the lump 
Third, Number of lumps detected. 
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Figure 4.27 Plot X,Y vs Z for a CGN at 45kft. 



POSITION 




a. Procedure to Detect the Lump 

1. Check a big or small lump by testing the conditions 
of three varying values of Cy increments (ACy) in 
sequence. If they are positive it is a big lump and 
set the flag- lump to 1, otherwise it is a small lump 
and set the flag- lump to zero; then go to 2 

2. Check the present knot position to see where its 
Z value equal to zero or to maximum Z value, if it is 
Z maximum then stop, it not go to 3 

3. If the process begins at the first knot position, set 
the begin and the flag-end to zero; check the status 
of the flag-lump for 1 (big lump) or 0 (small lump), 
if it is a big lump, then go to 4. If it is a small 
lump, then go to 7. 

4. Check the status for the begining or ending of the 
lump. If the flag-begin is 1, it represents that the 
beginning of a lump is found, then go to 5. Otherwise 
it's not found, then go to 6. 

5. Find the ending of the big lump by testing the 

conditions of 3 values of Cy increments in sequence. 

The first 2 should be negative and the third should 
be constant or positive. If the condition are 
satisfied, store the Z value of the position of the 
third knot, and the flag-begin to zero; then go to 
10 . 

6. Find the beginning of the big lump by testing the 

conditions of 3 different values of Cy increment 

(ACy) in sequence. The first Cy should be negative 
or constant, the second should be positive, and the 
third should be positive. If the conditions are 
satisfied, store the Z value of the position of the 
second knot, and set flag-begin to 1; then go to 10. 

7. Check the status for the beginning or ending of the 

lump. If the flag-begin is 1, it represents the 
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beginning of a lump is found, then go to 8. Otherwise 
it's not found, then go to 9. 

8. Find the ending of the small lump by testing the 

conditions of 3 values of Cy increment (A Cy) in 
sequence. The first one should be negative or 

constant. If the conditions are satisfied, store the 
Z value of the position of the second knot, and set 
the flag-begin to zero; then go to 10. 

9. Find the beginning of the small lump by testing the 
conditions of 3 values of Cy increment (A Cy) in 
sequence. The first should be constant or negative, 
the second should be positive, the third should be 
constant. If the conditions are satisfied, store the 
position of the second knot, and set the flag-begin 
to 1; then go to 10. 

10. Move to the next knot, then go to 2. 

This procedure is shown in Figure 4.29 and the detail of 
each procedure is shown in Appendix D. 
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CKECK AT THIS KNOT WHETHER 
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Figure 4.29 Flow Chart for Detecting Lumps. 



2 . To Determine the Area Under a Lump 

The area under the lump can be determined by 

AAREA = AC|ACy +CyACy (4.10) 

Suppose there is a lump as shown in Figure 4.30_. 
The area increment can be calculated by using eq(4.10). Then 
the area under BC(which is negative) is added to the area 
under AB. 




Cx 



Figure 4.30 First Procedure to Determine the Area. 

Next, area under CD is calculated, which is positive 
and add this to the last resulting area. Obviously, the sum 
would represent the total area of the lump as shown in 
Figure 4.31. 
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APQ-PBC 



Figure 4.31 Step 





CLCDE 



Step Procedure to 





Determine the Area. 



F. CONSTRUCTION OF THE DECISION TREE 

The characteristics of ships used for classification 
are, the number of lumps, the area of lumps, the knot posi- 
tions^ value) and where the beginning of the lump and lump 
maxima relative to the midships are located. 

In view of the above characteristic, it is necessary to 
find the relationships among them to make classification 
possible. Therefore, for each of the eight different class 
of ships, the decision tree is constructed which is based 
upon relationships observed in the plots of the knot posi- 
tion (Z values) for the beginning of the lumps normalized by 
the total ship length (Z value) VS. the number of lumps; the 
area of the lumps normalized by the total ship length 
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(2 value) squared vs the number of lumps; and the Z value of 
the lump maxima, as shown in Figure 4.32 through 
Figure 4.39. 

Thus, according to the number of lumps presented in the 
profile used as the first criteria, ships can be divided 
into 3 distinct groups. Then, classification, can be made 
by comparing further characteristic as shown in Table III, 
and the complete decision tree constructed from Table III is 
shown in Figure 4.40. 
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Figure 4.32 Plot the Beginning, Area, and Peak of a FF. 
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Figure 4.33 Plot the Beginning, Area, and Peak of a LST. 
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Figure 4.34 Plot the Beginning, Area, and Peak of a DD . 
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Figure 4.35 Plot the Beginning, Area, and Peak of a DDG. 
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Figure 4.36 Plot the Beginning, Area, and Peak of a 
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Figure 4.37 Plot the Beginning, Area, and Peak of a AOR. 



95 



BEGIN FROM MIDSHIP 




AREA OF LUMP 

O.OI 



0.005 



*3 

z 



S -50 , N -26 

3 = FREIGHTER AT RANGE 53000 FEET 



5 = 50 



V= 24 



3 NO. OF LUMP 



AREA OF LUMP 




PEAK FROM MIDSHIP 



Figure 4.38 Plot the Beginning, Area, and Peak of a Freighter. 
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Comparison of Different Types of Ships 
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Figure 4.40 The Decision Tree. 



G . SUMMARY 

The decision tree which is constructed from Table III 
does not provide 100-percent correct classification for all 
types of ship images. Furthermore, the presence of noise in 
images may cause complicatations in classifying the ships. 
For example, with excessive noise in the original image, 
Sobel operator for edge enhancement in the preprocessing 
process, still yield result with residual noise present in 
Figure 4.41. These residual noise is undesirable since it 
causes failure to extract profiles which retain necessary 
informations from the original images. The subsequent clas- 
sification of profile by the Fourier Coefficient method and 
the B-spline Coefficient method becomes difficult. 




Figure 4.41 Noisy Image. 

As discussed before in order to reduce the noise, an 
appropriate threshold gray value for the preprocessed image 
is set. This results in a silhouette image. However, if the 
value is too high the profile becomes broken which prevent 
successful operation in the closing process discussed in 
chapter 2, as shown in Figure 4.42. On the other hand 
decreasing the threshold value results in erroneous profile, 
as shown in Figure 4.43. In some cases, when the closing 
process is used, certain vital information is lost. This 
effect can be seen in comparing the image in Figure ‘‘.44, 
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and Figure 4.45. Therefore, loss of informations due to the 
attempt of eliminating noise and the inability in the 
preprocessing process to cope with the residual noise, yield 
erroneous profile. Consequently, failure occurs in the 
succeeding classification steps. 




Figure 4.42 High Threshold. 
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Figure 4.44 



Before Closing. 




Figure 4.45 After Closing. 
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V. CONCLUSION 



The shape of the superstructure of a ship is the most 
important feature used in the classification algorithms. To 
extract the edge profile, a Sobel operator is employed. The 
limitation of a Sobel operator is that some small details of 
the edge is lost. For example, the image of a DD at a range 
of 79000 feet after applying the Sobel operator shows the 
superstructure and the radar. But the edge of a small mast 
disappeared as shown in Figure 5.1. Furthermore, if the 
threshold value is set too big, the edge image of the super- 
structure profile becomes broken. The top superstructure 
profile is obtained by setting the gray value under the 
slope between the bow and the stern to zero. We need to 
apply a contour tracking process to refine the superstruc- 
ture profile. For some images, the connection of the broken 
profile pieces may be achieved in a Closing operation as 
discussed in Chapter 2. The disadvantage of the Closing 
operation is that some small details of the profile may 
disappear. The superstructure profile of a freighter at a 
range of 53000 feet is shown in Figure 5.2. After the 
Closing process the detail of the cranes disappeared as 
shown in Figure 5.3. 

The ship classification can be achieved by either the 
Fourier Coefficient method or the B-spline Coefficient 
method. In using these coefficients to classify ships, it 
is found that only the initial coefficients that lie between 
the 0th to the 20th, are relevant, while the rest are not. 
Inspection of the comparison curves of the same class of 
ships shows that, similarities in patterns exists up to the 
20-th point. Beyond that, diversities in shape are so great 
that inclusion of those additional points for classification 
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will be of no use. 
Figure 5.5. 



Examples are shown in Figure 5.4 and 




Figure 5.1 Edge Image of a DD at a Range of 79000 feet. 




Figure 5.2 Before Closing Process. 




Figure 5.3 After Closing Process. 
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Figure 5.4 Logarithmic Magnitude of a CGN at a Range of 45000 feet. 
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Figure 5.5 Logarithmic Magnitude of a CGN at a Range of 55000 feet. 



In the B-spline coefficient method the technique of 
uneven knot selection has been employed. With the original 
sampling points on the ship profile as input, a smaller set 
of approximation samples are collected. These can be used to 
reconstruct the ship profile with sufficient information 
retained from the original image. The execution time as 
approached to that of handling the original data directly 
has been greatly reduced. The reduction of the number of 
sample points by almost a factor of 10 is common. For a CGN 
ship a set of original sampling points of 290, has been 
reduced to 36. 

Comparing the two methods, Fourier Coefficient and 
B-spline Coefficient, the former method, in some cases is 
not effective to establish satisfactory classification of 
ships. This is due to difficulties in matching similarities 
of the shape of the coefficient curves. The latter, 
however, surpasses the former in that it is able to classify 
more ships accurately using computer programs. It is 
possible to improve the reliability of those two methods by 
reducing noise in the data collection process and the 
preprocessing process. 
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APPENDIX A 



THE PROGRAM TO OBTAIN THE SUPERSTRUCTURE 
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6-Dec-1984 17:13:40 
Source Listing 6-Dec-l98* 17:18:31 

program cut(input,output»infile,outfilo) ; 

type 

byte = 0..255; 

lmagarotul = packed array CO. .257 3 of byte? 
imagaroui2 - packed array CO. . 255 3 of integer; 
imagerow3 =■ packed array CO*. 2 55 3 of byte; 
roml = oacked array CQ-.2553 of byte; 
ro«2 s oacked array CO*. 2553 of integer; 

var 

sobel : array CO. .633 of iraa gerow2 ; 

l * J : byte; 

f : array CO. .653 of imageroilj 
outfile : file of imagero<u3; 

d* * dy t range * br ight t n , m Jinteger; 
infile : file of roiul! 

slope : real ; 

image larray CO. . 633 of rouii; 
x t *1 »x2 t y 1 • y2 : integer; 

NUM1, NUM2,NUM3,NUM4 t max,min # (»axl ; integer; 
cal :array CO. . 63 3 of rouii ! 
thes : integer ; 

NAME : PACKED ARRAT Cl. . 20 3 OF CHAR; 

BEGIN 

WRITELNC ‘INPUT SHIP FILENAME OUT CUT. DAT*); 

RE ADLN(NAME) ; 

WRITELNC •THRESHOLD*) ; 

RE ADCTHES) ; 

WRITELNC * NUM TO CUT FROM LcFT*); 

RE AO(NUMl) : 

WRITELNC* N JM2 TO CUT FROM RIGHT*); 

REA0CNUM2) ; 

WRITELNC *NUM3 TO CUT FROM TOP*); 

RE A0CNUM3) ; 

WRITELNC * NUM4 TO CUT FROM BUTT0H*); 

R EAO CNUM4 ) ; 

open C in f i 1 e * NAME » hi stor y : = old t 
access.method .^sequential i 

record.length : = 256» record_type : = f ixed) ; 
open C ou t f i 1 e « * C UT . da t * » hi s t or y : =neui » re cord _ 1 eng t h :=256 » 
record_type *. = f ixed) 5 
resetCinfile) ; 
rewrite (outfile); 
i : = 0 ; 

while not eof Cinfile) do 
oegin 

read C in f i 1 e , i ma ge C i 3 ) l 
for j : = 0 to 255 do 

fCi*l*j+13 := imageCitj3; 
i 2 *i ♦! ; 
end; 

(^compute sobel$) 
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Source Li s ting 



6-0ec-196<* 17:18:40 

6-Dec-1984 17:18:31 



V A X- 
_UR t 



FOR J : =0 TO 7 00 

c C64*J*2*S3:=FC63*J*2483: 
for i := 0 to 63 do begin 
f C i ♦ 1 . 0 3 := f C i ♦ 1 • 1 3 ; 
fCi*l,2573 :»fCi*l .2563; 
end ; 

for j:= 0 to 255 do begin 
fCO.j+13 := fCl.j+13: 
fC65,j*13 := fC64.j+13; 
end ; 

f C3. 03 : = f C 1, 13 ; 
fC0.2573: = fCl.25631 
fC65*03:=fC64 t 13; 
f C65.2573:=fC64,2563 ; 

($set initial max. min $3 

ma x : =0 ; 
maxi • — 0 * 
rain : =25 5 : 

for i:= 0 to 63 do begin 
for j:= 0 to 255 do begin 

dx:=fCi,j3*2-fCi*l.j3+fCi*2.j3-fCi.j*23 
-2*f Ci*l • j*23-f Ci*2 • j + 2 3; 
dy: = fCi*2,j3 + 2*fCi*2.j-U3 + fCi-*2,j + 2 3-fCi.j3 
-2*fCi» j + D-f Ci. j + 23: 
sobelCi»j3:=round((dx*S2+dy*$2)*30.5); 
if max < sobelCi.jJ then 
max := sobe 1 Ci , j 3 ; 
if rain > sobelLi»j3 then 
min := sobelCi*j3; 

end ; 
end : 

range J= max-min? 

($ rescalot) 

for i : = 0 to 63 do begin 
for j:= 0 to 255 do begin 

calCi.jO := r oun d ( ( C sob e 1 C i « j 3-m in )£ 2 55 ) / r ange 1 i 
if (j<-NUMl ) or ( J>=255-NUM2) then 
calCi»j3 : =0 ; 

if (i<=NUM3) or (i>=63-NUKO then 
ca 1C i • j 3 := 0: 

($ cut thesholdt ) 

IF TMESOO THEN BEGIN 
if calCi»J3 <-thes 
then 

calCi.j3 0 
else 

c a 1 C l » j 3 := 255 ; 

end; 

end; 
end ; 

(Sfind point at raw and aft$) 



no 



Source Listing 



6-0ec-1984 17:13:40 

6-Dec-l 964 17:18:21 



in : = o; 
n : = o : 
bnght:= 255: 

for j := 0 to 150 do begin 
for i : = 0 to 63 do begin 

l f < c a 1 U l * j 3 = br i gh t ) and (n=0) then 
begin 

y 1 : = j ; 

x l : - i : 
n n+ l : 
end; ( $i f ♦ ) 

i f CcalCi * 255- j3 = bright ) and (m = 0) then 
begin 

y 2 := 2 55-j; 

x 2 : = i ; 

v) : = «♦ i ; 

end; ( ♦ l f ♦ ) 

end;(-forv) 

endlC^fori) 

»riteln('xl = 'ixl, "yl='*yl» "x2 = '*x2»'y2=‘’’»y2) ; 

(* cut line« ) 

slooe : = 1 . 0 : 
if x 1 » x 2 then 
begin 

slooe : = 0 ; 

for i:=xl»l to 63 do begin 
for j : = y 1 - 5 to y2«-5 do 
c al C i • j 3 : = 0: 
end; ( ♦ f or* ) 
end ; ( 3 i f ♦ ) 
i f slopeOO then 
begin 

slope :=(x2-xl )/( y2-y 1 ) : 
if slope < 0 then begin 

slope := absCslope); 
for j := y 1 to yZ do begin 

x := x 1 ♦ 1-round ( ( j-yl )<= slope ) ; 
for i :=x to 63 do 
calCi*j3:=o; 
end : ( £ f o r ♦ 3 
end; (Sif$) 

if slope >0 then begin 

for j:* yl to y2 do begin 

x :* x 1 ♦! ♦round C ( j-y 1 )«slope) ; 
for i : = x to 63 do 
calCi » j3 :=0 : 
end: (*for?) 
end ; (Si f ♦ ) 
mri teln( # step2 * ) ; 

•nd; (iif*) 

(♦put ou t f i 1 et ) 

for i:=Q to 63 do begin 
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ill 



Sourc e Listing 



6-D9C-1984 17:13:40 

6-Dec-1984 17:18:31 



VA X- 



for j : = 0 to 255 do 

outfiie A Cj3 := calCi,j3: 
put Coutfile): 
end : C^for*) 
closeCin f i le) : 
end . 
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APPENDIX B 



THE PROGRAM OF CONTOUR FOLLOWING 
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10-Dec-1984 13:53:33 
Source Listing 6-0ec-l984 17:33:13 

PROGRAM TR AC ECINPUT. OUTPUT. INFILE. OUTFILE): 

TYPE 

BYTE = 0..255; 

IMAGEROW1 = PACKED ARRAY CO. .2 55 3 OF BYTE: 

RDWI = PACKED ARRAY CO. . 2573 OF BYTE; 

VAR 

R*C*F0*Fl.F2.F3*F4,F5»F6.F7.F8 : BYTE: 

F : ARRAY CO. .653 DF RDWi: 

A, IMAGE : ARRAY CO. .63 3 DF IMAGE RDWI ; 

infile : file of imagerdui: 

R1 ,C1 , R2.C2 .COIR, I , J. COUNT : INTEGER: 

ROW, COL : ARRAY CO.. 5123 OF INTEGER: 

dutfile : file df imagerowi: 

NAME : PACKED ARRAY Cl. .203 OF CHAR; 

procedure store; 



BEGIN 

CDLCI 3 : -C- 1 : 

R0wci3:=R-i: 

WRITELNC*CDL=' •COLCI3, *ROW=* f ROWCI3): 
count : = i : 
l:=I*i: 
eno: 



PROCEDURE CMOVE; 



BEGIN 

FO : = F C R- 1 tC 3 : Fi: = FCP-l,C*13; F2:=FCR.C+13: F 3 : = FCR*1 .C*l 3 ; 
F4:«FCR^1.C3! F5:=FCR^1*C-13: F6:*FCR,C-13: F7:«FCR-1.C-13: 
CASE COIR OF 
O: BEGIN 

*. ENO 



IF C F 3 = 


0 ) AND 


C F 2 = 


255) THEN 


BEGIN 


C : =C 


♦i: COIR 


ELSE IF 


C F 2 = 0 ) 


ANO 


(F1=2S5 ) 


THEN 


BEGIN 


R : =R- 1 ; 


STORE 


: c : = c ♦ 1 ; c di r : = 0 : eno 






ELSE IF 


( F 1 = 0 ) 


AND 


CF0=255) 


THEN 


BEGIN 


r : = r - 1 : 


COIR : 


=0: ENO 












ELSE IF 


C F 0 = 0 ) 


AND 


CF7=255) 


THEN 


BEGIN 


0 

11 

0 

1 


STORE 


: r:=p- 


1 : coir:=3: eno 






ELSE IF 


C F 7 = 0 3 


ANO 


C F 6 = 2 5 5 ) 


THEN 


BEGIN 


c : = c-i : 


COIR: 


=3: END 












ELSE IF 


( F 6 = C ) 


AND 


( F 5=2 5 5 ) 


THEN 


BEGIN 


R : = R ♦ 1 ; 



STORE: C :=C“1 : COIR:=3; END 

else begin r:=r*i: coi r : = 2 : eno: 
eno: 

l: BEGIN 

IF ( F 5 = 0 ) ANO (F4= 255) THEN 6EGIN R:=R*11 



C0IR:=2: END 

ELSE IF (F4=0) and CF3=255) THEN 

store; r:=R*i; coir:=i; eno 

ELSE IF ( F 3 = 0 ) AND (F2 = 255 ) THEN 
CDIR : = 1 : ENO 

ELSE IF ( F 2 = 0 ) AND (Fl=255) THEN 

store: c:=c*i: coir:=o: end 

ELSE IF ( F 1=0) AND (F0 = 255 ) THEN 

cdir:*o; eno 



BEGIN C : = C ♦ 1 ; 
BEGIN C:=C*i: 
BEGIN R : = R - 1 ; 

begin r:=r-i: 



VAX- 11 
_ QR A 0 : 
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iO-Dec-1984 13:53:33 
Source Listing 6-Dec-198*t 17:38:13 

ELSE IF C F 0 = 0 ) AND (F7=255> THEN BEGIN C:=C-i: 

store; r : = r - i ; cdir:=o; end 
else BEGIN C:=c-i; CDIR : =3 ; end; 
end ; 

2: BEGIN 

IF ( F 7 = 0 ) AND ( F6 = 2 5 5 ) THEN BEGIN C 2 *C -1 I 
CDIR:=3; END 

ELSE IF ( F 6 = 0 ) AND (F5=255) THEN BEGIN R:=R*i; 

store; c:=c-i: cdir:=3; end 

ELSE IF C F 3 = 0 ) AND CF4=255; THEN 6 EG I N R:=R*i; 

C D I R : = 2 ; END 

ELSE IF ( F 4 - 0 ) AND (F3=255) THEN BEGIN C : = C ♦ 1 2 

store: r : =r> i ; cdir:=i; end 

ELSE IF ( F 3 = 0 ) AND (F2 = 255 ) THEN BEGIN C:=C*15 
C D I R : = 1 ; END 

ELSE IF ( F 2 = 0 ) AND (Fl = 255 > THEN BEGIN P:=R-i; 

store; c:=c*i; cdir;=i; end 
else begin r:=r-i; cdir:=0; end; 
end; 

3: BEGIN 

IF ( F 1 = 0 ) AND CF0 = 255) THEN BEGIN R : = R - 1 ; 

cdir:=o: end 

ELSE IF ( F 0 = 0 ) AND (f=7 = 255 ) THEN BEGIN C : - G — 1 5 

store: r : = p-i ; cdir:=o; end 

ELSE IF ( F 7 = 0 ) AND (F6 = 255 ) THEN BEGIN C:=C-i; 

CDIR: =3; END 

ELSE IF ( F 6 = 0 ) AND (F5 = 255 ) THEN BEGIN R:=R*i: 

store; c:=c-i; cdir:=2; end 

ELSE if ( F 5 = 0 ) AND (F4=255) THEN BEGIN R:=R*i; 

C DI R : = 2 : END 

ELSE IF CF4=0) ANC (F3=255) THEN BEGIN C:=C*i: 

store; r:=r*i; cdir:=2; end 

ELSE BEGIN C : = C + 1 ; CDIR:=i; END; 

enj ; 
end; 
eno; 



PROCEDURE initial; 

BEGIN 

FOR I ;=0 TD 255 DO BEGIN 

cqlcid :=o; 

R owe 1 3 :=o; 

end; 

FOR J : *0 TD 63 DD BEGIN 
FDR I : =0 TO 255 DD 
AER,C3:=0; 

end; 

I : = 0 ; 
end; 

PROCEDURE FIRST; 

VAR 

n»m :integer; 

BEGIN 



VAX’ 11 
_ UR A o : 
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Source Lasting 



10-Oec-19bL 13:53:33 
6-Dec-198L 17:38: 13 



N : = 0 : m:=o; 

FOR C : =0 TO 200 00 BEGIN 
FOR R : = 0 TO 63 00 BEGIN 

IF ( FC R , C 3 =25 5 ) ANO (N = 0) THEN BEGIN 
Cl : = c; R 1 : = R ; n : = n+i; 

ENO; 

IF (FCR,255-C3=255) ANO (M = 0) THEN BEGIN 
C2 := 255-C: R2 : = R; m :=m*i; 

eno: 

eno: 

Eno: 

WRITE LN('COLl=*.Cl,'ROWl=',Rl t 'COL2=',C2,*ROW2='.R2); 

eno; 

(*$* *2* *<:****$<:-<£ ***;5:*$SS ) 

(* MAIN PROGRAM $) 

BEGIN 

WRITELNC 'INPUT CUT OR CONLINE FILENAME') : 

R EAOLNC NAME ) ; 

OPEN ClNFlLEtNAME.HlSTORY := OLO, 

ACCESS. METHOD : = S E 3UENT I AL • 

RECORO. LENGTH : = 256 • RECORO.T YPE :=FIXEO): 

OPEN CDUTFILE, 'TR.OAT'.HISTORY : = N E U , R c C □ R 0_L E NG TH :=256. 
RECORO.T YPE : =FI X EO ) ; 

resetcinfile) : 

REWRITE (QUTr ILE) ; 

R : = 0 : 

WHILE NOT EOF CINFILE) 00 
BEGIN 

READ (INFILE, IMAGEER3) ; 

FOR C : - 0 TO 255 00 

FCR*l,cn3 : = IMAGED RiC3: » 

R := R ♦ 1 ; 

END; 

first; 

r :*ri: c :=ci: coir :-i: 
initial; 

while CCOC2) 00 BEGIN 
store; 
cmoye ; 
eno; 

FOR i:=0 TO COUNT 00 

AC ROWE I3»C3LCI33 1=255 ; 

FOR R : =0 TO 63 00 BEGIN 
FOR C:=C TO 255 00 

OUTFILE*CC 3 := A t R * C 3 I 
PUT < OUT F I L E ) : 

eno; 

WRITELNC'NUMQ6R=', COUNT); 

closecin c ile); 

ENO. 



VAX-11 
_QR AO : 
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APPENDIX 
THE SUBROUTINE 



C 

"PARAM" 
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SUBROUTINE PARAM(X,Y»Z,W f KfZB,ZE»K,S,N,T,CX,CY,FP,IOPT,IPAR,IER) 0930 

C GIVEN THE SET OF DATA POINTS (X(I),Y(I)) WITH CORRESPONDING Z- 0940 

C VALUES ZCI). 1=1*2, M AND GIVEN ALSO THE SET OF POSITIVE 0950 

C NUMBERS W(I),I=1,2, ...M, SUBROUTINE PARAM FINOS A SMOOTH APPROXIMAT- 0960 

C ING CURVE WITH PARAMETER REPRESENTATION X = SX(Z), Y = SY(Z). 0970 

C $X(Z) ANO SY(Z) ARE TWO SPLINE FUNCTIONS OF DEGREE X WITH THE NUMBER 0980 

C AND THE POSITION OF THE KNOTS T ( J ) , J = 1 , 2 * . . . N AUTOMATICALLY 0990 

C CHOSEN 3 Y THE ROUTINE. THE SMOOTHNESS OF SX(Z) AND SY(Z) IS 1000 

C ACHIEVED BY MINIHALIZING THE SU M( 0 X ( R )**2 *0 Y ( R ) **2 > WHERE DX(R) 1010 

C ANO 0 Y ( R ) STAND FOR THE 0 I SC ONT I NU I T Y J U M P OF THE K T H DERIVATIVE 1020 

C OF SX(Z) ANO S Y ( Z ) AT ThE KNOT T ( R ) f R=K*2, N-K-l. 1030 

C THE AMOUNT OF SMOOTHNESS IS DETERMINED BY THE CONOITION THAT F(P) * 1040 

C SUM(W(I)*((X(I)-SX(Z(I)))**2 ♦CYCX)-SYCZCI)))«*Z)) BE <= S. WITH 1050 

C S A GIVEN NON-NEGATIVE CONSTANT. 1060 

C THE SPLINE FUNCTIONS SX(Z) ANO SY(Z) ARE GIVEN IN ThEIR B-SPLINE 1070 

C REPRESENTATION (B-SPLINE COEFFICIENTS CX(J) t RESP. CY(J),J = 1, N-K-l) 1030 

C ANO can BE EVALUATEO BY MEANS OF FUNCTION OERIV. 1090 

C CALLING SEQUENCE: 1100 

C CALL PARAMCX, Y,Z,U,M, Z3,ZE,K f S , N , T , C X , C Y , F P , ICPT, IPAR* IFR) 1110 

C 1120 

C INPUT PARAMETERS: 1130 

C X : ARRAY, LENGTH M, CONTAINING THE ABSCISSAE of THE OATA POINTS 1140 

C Y : AR R A Y , L ENGTH M, CONTAINING THE ORDINATES OF ThE DATA POINTS 1150 

C W : ARRAY, MINIMUM LENGTH M, CONTAINING THE WEIGHTS W(I). 1160 

C M : INTEGER VALUE, CONTAINING THE NUMBER of OATA POINTS. 1170 

C K : INTEGER V A LU E , C 0 NT A I N I NG THE OEGREE OF SX(Z) ANO SY(Z). 1180 

C S : REAL VALUE, CONTAINING THE SMOOTHING FACTOR. 1190 

C IOPT : INTEGER VALUE WHICH TAKES THE VALUE 0 OR 1. 1200 

C 1 0 P T = C : THE ROUTINE WILL RESTART ALL COMPUTATIONS. 1210 

C IOPT=i: THE ROUTINE WILL START WITH THE KNOTS FOUND AT THE 1220 

C LAST CALL OF THE ROUTINE. IF I0PT=1 THE OUTPUT 1230 

C PARAMETERS T ANO N ARE INPUT PARAMETERS AS WELL. 1240 

C IF 1 0 P T = 1 THE USER MUST PROVIOS WITH A COMMON BLOCK 1250 

C COMMON/OPT1/NROATA(NEST),FPO,FPOLO,NPLUS 1260 

C IPAR : INTEGER flag. 1270 

C IPAR = 0: FOR EACH OATA POINT (X(I),Y(I)) THE PROGRAM AUTOMATICALLY 1280 

C CHOOSES A CORRESPONDING VALUE OF THE PARAMETER Z, I.E. 1290 

C ZCl)*0:ZCI)>ZCI-l>*SQRT(CXCI)-XCI-l)>s*2«CYCI)-YCI-l))*«2) 1300 

C THE BOUNDARIES FOR THE PARAMETER Z ARE CHOSEN AS FOLLOWS 1310 

C ZB =■ 2(1) : ZE = Z ( M ) . 1320 

C IPAR = l: THE USER HIMSElF PROVlOES WITH THE VALUES OF ThS 1330 

C PARAMETER Z ANO WITH THE BOUNOARIES ZB ANO ZE. 1340 

C Z : array, LENGTH M, CONTAINING THE VALUES OF THE PARAMETER Z 1350 

C (IPAR = 1) 1360 

C Z8,ZE: REAL VALUES, CONTAINING THE BOUNOARIES OF the PARAMETER Z 1370 

C (IPAR =• 1). 1380 

C 1390 

C OUTPUT PARAMETERS: 1400 

C T : ARRAY, LENGTH NEST (SEE OATA INITIALIZATION STATEMENT), 1410 

C WHICH CONTAINS THE POSITION OF THE KNOTS, I.E. THE POSITION 1420 

C OF THE INTERIOR KNOTS TCK+2), T(N-K-l), AS WELL AS THE 1430 

C POSITION OF THE KNOTS T ( 1 ) = T ( 2 )=...= T ( K ♦ 1 ) = Z B AND ZE = 1440 

C T(n-K)=...sT(N) WHICH ARE NEEOEO FOR THE B-SPLINE REPRESENT. 1450 

C CX,CY: ARRAYS, LENGTH NEST, CONTAINING the B-SPLINE COEFFICIENTS 1460 

C OF S X ( Z ) , RESP. SY(Z). 1470 

C N : INTEGER V A LU E , CO N T A I N I NG THE TOTAL NUMBER OF KNOTS. 1480 

C FP : REAL VALUE, WHICH CONTAINS THE S UM ( W I * ( x I - S X ( Z I ) ) **2 ) 1490 

C ♦ SUM(WI*CYI-SY(ZI))**2), 1=1, 2, ...M. 1500 

C IER : ERROR COOE 1510 

C I E R = 0 : NORMAL RETURN. 1520 

C I E P =- 1 I NORMAL RETURN f SX(Z) AND SY(Z) ARE I N T F R P JL A T I N G SPLINES 1530 
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oo<”>noonooon<”>or»o onoooooooooo ooooooooor>oonooor>onr»no 



IER=-2:nORHAL RETURN, SX(Z) ANO SY(Z) are POLYNOMIALS OF OEGREE K 1540 

I E R > 0 .-ABNORMAL TERMINATION- 1550 

IER= 1 : THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE 1560 

STORAGE SPACE, SPECIFIED BY THE PARAMETER NEST. 1570 

PROBABLY CAUSESINEST OR S TOO SMALL. 1530 

IER=2: A THEORETICALLY IMPOSSIBLE RESULT WAS FOUND DURING 1590 

THE ITERATION PROCESS. 1600 

PROBABLY CAUSESITOL TOO SMALL 1610 

IER=3:THE MAXIMAL NUMBER OF ITERATIONS MAXIT HAS BEEN 1620 

REACHED. 1630 

PROBABLY CAUSESiMAXIT OR TOL TOO SMALL. 1640 

IERMOISOME OF THE INPUT DATA ARE INVALIOCSEE RESTRICTIONS). 1650 

1 660 

RESTRICTIONS: 1670 

1) M > K > 0 1680 

2) 2B <= Z ( R ) < Z ( R ♦ 1 ) <= ZE, R = 1 , 2 , . . . M- 1 . (IPAR * 1) 1690 

3) W(R) > o, R= 1 , 2 , . . . M. 1700 

4) S >= 0. 1 MO 

5) NEST >= 2*X42. 1720 

1730 

OTHER SUBROUTINES REQUIRED: 1740 

eSPL I N,COSSIN, ROTATE, BACK, NK NOT, DISCO ANO RATION. 1750 

DIMENSION X(M),YCM),WCM),2(M),T(200),CX(2G0).CYC200), 

< FPINT(200),RX(200),RY(200),OIAG(200),OPRIMEC200), 

< GC20C,6),B(200*7),QC400,6),H(7),NR0ATA(200),A(200.5) 

C0MMCN/0PT1/NRCATACNEST ) ,FPO,FPOLO,NPLUS 1790 

NROATA: INTEGER ARRAY, LENGTH NEST, WHICH GIVES THE NUMBER OF 1800 

DATA POINTS INSIDE EACH KNOT INTERVAL. 1810 

FPO : REAL VALUE, WHICH CONTAINS ThE S UM C W I s( X I- S X C Z I > ) s*2 ) ♦ 1820 

SUM( WI*C YI-SYCZI >)**2) WITH SXCZ) ANO SY(Z) L E A S T - S 3U A R E S 1830 

POLYNOMIALS OF OEGREE K. 1840 

F POL 0 : REAL VALUE, WHICH CONTAINS THE SUM ( W I * C X I- S X ( Z I ) ) ** 2 ) ♦ 1650 

SUMCWI*( YI-SYCZI))**2) WITH SX(Z) AND SY(Z) LEAST- SOUAR E S I860 

SPLINE FUNCTIONS CORRESPONDING TO THE LAST FOUNO SET OF 1870 

KNOTS BUT ONE. 1880 

NPLUS : INTEGER VALUE, GIVING THE NUMBER OF KNOTS OF THE LAST 1B90 

SET MINUS THE NUMBER OF THE LAST SET BUT ONE. 1900 

C0MM0N/0PT1/NR0ATA, FPO, FP OLD, NPLUS 1910 

DATA INITIALIZATION STATEMENT TO SPECIFY 1920 

TOL : THE REQUESTEO RELATIVE ACCURACY FOR THE ROOT O c FCP) = S. 1930 

MAXIT: THE MAXIMAL NUMBER OF ITERATIONS ALLOWEO. 1940 

NEST J AN OVER-ESTIMATE OF THE NUMBER OF KNOTS N. THIS PARAMETER 1950 

MUST BE SET BY THE USER TO INOICATE THE STORAGE SPACE 1960 

AVAILABLE TO THE SUBROUTINE. THE OlMENSION SPECIFICATIONS 1970 

OF THE ARRAYS T , C X , C Y , N R 0 A T A , F P I N T , R X , R Y , 0 1 A G , 0 P R I H E ( N ) , 1930 

A(N,K),G(N,K*1),B(N,K*2),2(M,K*1) ANO H(K*2) OEPENO 1990 

ON N , M ANO K. SINCE N IS UNKNOWN AT THE TIME THE 2000 

USER SETS UP THE OlMENSION INFORMATION AN OVER-ESTIMATE 2010 

OF THESE ARRAYS WILL GENERALLY BE MAOc. THE FOLLOWING 2020 

REMARKS ARE INTENOEO TO HELP THE USER 2030 

1) 2*K+2 <= N <= M+K+I 2040 

2) THE SMALLER THE VALUE OF S, THE GREATER N WILL BE. 2050 

3) NORMALLY N - M/2 IS AN OVER-ESTIMATE. 2060 

DATA TOL/O. 001/, MAXI T/20/, NEST/200/ 

C BEFORE STARTING COMPUTATIONS A DATA CHECK IS MAOE. IF THE INPUT 2080 

C DATA ARE INVALIO CONTROLE IS IMMEDIATELY REPASSEO TO THE ORIVER 2090 

C PROGRAM CIEP-10). 2100 

IER = 0 2110 

K 1 = K ♦ 1 2120 

K 2 * KIM 2130 

NMIN = 2*K1 7 1 4 C 
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IFCK.LE.3) IER = 10 2150 

IFCM. LT.K1 .OR. NEST.LT.NMIN) IER - 10 2160 

IFCS.LT.O.) IER = 10 2170 

IFCIER-NE.O) GO TO 440 2190 

C CHECK WHETHER THE Z-VALUES ARE PROVIOEO WITH BY THE USER. 2190 

IFC IPAR.NE.O ) GC TO 6 2200 

C FINO FOR EACH DATA POINT A CORRESPONDING VALUE OF THE PARAMETER Z 2210 

C AND FIX THE BOUNDARIES ZB ANO ZE. 2220 

ZCI) = 0. 2230 

00 4 I = 2 v M 2 240 

ZCI) = ZCI-1H'SQRTCCxCI)-XCI-1))S*24CYCI>-YCI-1)>s*2) 2250 

4 CONTINUE 2260 

ZB = ZCI) 2270 

ZE * ZCH) 2280 

6 IFCZB.GT.ZC1) •OR. ZE.LT.ZCM) .OR. WCD.LE.O.) IER = 10 2290 

00 10 1 = 2. M 2300 

IF(ZCI-l).GE.ZCI) .OR. WCD.LE.O.) IER = 1 0 2310 

10 CONTINUE 2320 

IFCIER.NE.O) GO TO 440 2330 

C CALCULATION OF ACC, THE ABSOLUTE TOLERANCE FOR THE ROOT OF FCP)=S. 2340 

ACC = TOL*S . 2350 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2360 

PART l: DETERMINATION OF THE NUMBER OF KNOTS ANO THEIR POSITION C 2370 

^ 3 C 2390 

GIVEN A SET OF KNOTS WE COMPUTE THE LEAST- SO UAR E S SPLINES SXInFCZ) C 23 90 

ANO STINFCZI.IF THE SUM FCP=INF)<=S WE ACCEPT THE CHOICE OF KNOTS. C 2400 

OTHERWISE WE HAVE TO INCREASE THEIR NUMBER. C 2410 

THE INITIAL CHOICE OF KNOTS OEPENDS ON THE VALUE OF S AND IOPT. C 2420 

IF S = 0 WE HAVE SPLINE INTERPOLATION; IN THAT CASE THE NUMBER OF C 24 30 

KNOTS EQUALS NMAX = M*K*1. C 2440 

IF S > 0 ANO C 2450 

IOPT = 0 WE FIRST COMPUTE THE L E A ST- S QU A RE S POLYNOMIALS OF C 2460 

DEGREE k: N = NMIN = 2 s K ♦ 2 C 2470 

1 0 P T = 1 WE START WITH THE SET OF KNOTS FOUNO AT THE LAST C 2480 

CALL OF THE ROUTINE, EXCEPT FOR THE CASE THAT S > FPO; THEN C 2490 

WE COMPUTE DIRECTLY THE L E A S T- S QUA R E S POLYNOMIALS OF DEGREE K. C 2500 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2510 

C DETERMINE NMAX, THE NUMBER OF KNOTS FOR SPLINE INTERPOLATION. 2520 

NMAX = M ♦ K 1 2530 

IFCS.GT.O.) GO TO 45 2540 

C IF S=0, SXCZ) AND SYCZ) ARE INTERPOLATING SPLINES. 2550 

N = NMAX 2560 

C TEST WHETHER THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE ONE. 2570 

IFCN. GT.NEST) GO TO *20 2580 

C FINO THE POSITION OF THE INTERIOR KNOTS IN CASE OF INTERPOLATION. 2590 

MK1 = M-Kl 2600 

IFCMK1.EQ.0) GO TO 60 2610 

K 3 = K/2 2620 

1 = K 2 2630 

J = K3+2 2640 

I F C K 3s 2 . EC-K ) GO TO 30 2650 

00 20 L=1,MK1 2660 

TCI) = ZCJ) 2670 

I = 1 4 1 2680 

J = J 4 1 2690 

20 CONTINUE 2700 

GO TO 60 2710 

30 DO 40 L=1,MK1 2720 

TCI) = CZCJ )-ZCJ-l))*0.5 27 30 

I = I 4 l 2740 

J = J4l 
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40 CONTINUE 2 7 SO 

GO TO 60 2770 

C IF S>0 OUR INITIAL CHOICE OF KNOTS DEPEN03 ON THE VALUE OF IOPT. 2700 

C IF IOPT = 0 GR 1 0 P T = 1 ANO S>=FP0, WE START COMPUTING THE L E A S T - SO U A R E S 2 7 90 

C POLYNOMIALS O p DEGREE K WHICH ARE SPLINES WITHOUT INTERIOR KNOTS. 2800 

C IF I 0 PT = 1 ANO F P 0 > S WE START COMPUTING THE LEAST- SQUARE S SPLINES 2910 

C ACCORDING TO THE SET OF KNOTS FOUNO AT THE LAST CALL OF THE ROUTINE. 2820 

45 IFCIOPT.LE.O) GO TO 50 2e30 

IFCFPO.GT.S) GO TO 60 2840 

50 N = NMIN 2850 

NROATA(l) = M-2 2860 

C MAIN LOOP FOR THE OIFFERENT SETS OF KNOTS. M IS A SAVE UPPER BOUNO 2870 

C FOR THE NUMBER OF TRIALS. 2800 

60 00 200 ITER = 1 * M 2 890 

IF(N.EO.NMIN) IER = -2 2900 

C FINO NRINT, TNE NUM3ER OF KNOT INTERVALS. 2910 

NRINT = N— N M I N ♦ 1 2920 

C FIND THE POSITION OF THE ADDITIONAL KNOTS WHICH ARE NEEDED FOR 2930 

C THE B-SPLINE REPRESENTATION OF SXCZ) AND ST(Z). 2940 

NKI * N-Kl 2950 

I = N 2960 

00 70 J= 1 , K 1 29 70 

TCJ) = ZB 2980 

TCI) = ZE 2990 

I = I- 1 3000 

70 CONTINUE 3010 

C COMPUTE THE B-SPLINE COEFFICIENTS OF THE L E A S T- SQU A R E S SPLINES SXINFCZ) 3020 

C AND SYINFCZ). THE OBSERVATION MATRIX A IS BUILT UP ROW 3T ROW AND 3030 

C REOUCEO TO UPPER TRIANGULAR FORM 3T GIVENS TRANSFORMATIONS 3040 

C WITHOUT S3UARE ROOTS. AT THE SAME TIME FP=FCP=INF) IS COMPUTEO 3050 

FP = 0. 3060 

C INITIALIZE THE OBSERVATION MATRIX A. 3070 

00 80 1=1, NKI 3090 

DIAGCI) = 0. 3090 

RXCI) = 0. 3100 

RYCI ) = 0. 3110 

DO 80 J=1,K 3120 

AC I. J) » 0. 3130 

80 CONTINUE 3140 

L = K 1 3150 

00 1 30 I T = 1 , M 3160 

C FETCH THE CURRENT OATA POINT X C I T ) , Y C I T ) * Z C I T ) . 3170 

XI = XCIT) 3180 

Y I = YCIT) 3190 

ZI * ZCIT) 3200 

WI = WCIT) 3210 

C SEARCH FOR KNOT INTERVAL T C L ) <= ZI <= TCL*1 ). 3220 

IFCZI.GE.TCL+1) .ANO. L.NE.NK1) L = L*1 3230 

C EVALUATE THE C K ♦ 1 ) NON-ZERO 8-SPLINES AT ZI ANO STORE THEM IN Q. 3240 

CALL BSPLINCT,N, K,ZI ,L,H) 3250 

00 90 1 = 1, K1 3260 

IFCHCD.LT. 0.1E-07) H C I ) = 0. 3270 

OCIT, I) = HCI) 3230 

90 CONTINUE 3290 

C ROTATE THE NEW ROW OF THE OBSERVATION MATRIX INTO TRIANGLE BY 3300 

C GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS. 3310 

J = L-Kl 3320 

00 110 1 = 1, K1 3330 

IFCWI.EQ.o.) GO TO 130 3340 

J = J*1 3350 

PIV * HCI) 3 3t> 0 
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IFCPI V . E3. 0- ) GO TO 110 3370 

C CALCULATE the parameters OF THE GIVENS TRANSFORMATION. 3 3 B 0 

CALL COSSINCPIV, WI f OlAG(J) .COS. SIN) 3390 

C TRANSFORMATIONS TO RIGHT HANO SIOES. 3400 

CALL ROTATE(PIV,COS*SIN # XI*RX(J)) 34 10 

CALL ROTATE(PIV,COS*SIN*YI * R Y ( J ) ) 3420 

IFC1.EQ.K1) GO TO 120 3430 

12 = 0 3440 

13 = 1*1 3450 

DO 100 II = I 3 « K 1 3460 

12 = 1 2 + 1 3470 

C TRANSFORMATIONS TO LEFT HANO SIDE. 3430 

CALL RQTATECPIV,C0S.SIN.H<I1).ACJ.I2)) 3490 

100 CONTINUE 3500 

110 CONTINUE 3510 

C A 00 CONTRIBUTION OF THIS ROW TO THE SUM OF SQUARES OF RESIDUAL 3520 

C RIGHT HANO SIDES. 3530 

120 FP = FP*UI*<XI**2*YI**2) 3540 

130 CONTINUE 3550 

IFCIER.EQ.-2) FPO = FP 3560 

C BACKWARD SUBSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS. 3570 

CALL 3ACKCA*RX,NK1 « K v C X ) 3530 

CALL 3ACKC A f RY, NK1.K.CY) 3590 

C TEST WHETHER THE APPROXIMATION X = SX I NF < l ) , Y = S T I NF ( Z ) IS AN 3600 

C ACCEPTABLE SOLUTION. 3610 

FPMS = FP-S 3620 

IFCABS(FPMS).LT.ACC) GO TO 440 3630 

C IF F(P=INF) < S ACCEPT THE CHOICE OF KNOTS. 3640 

IFCFPMS.LT. 0.) GO TO 250 3650 

C IF N=NMAX,SXINF(Z) ANO SYlNFCZ) ARE INTERPOLATING SPLINES. 3660 

IFCN.EQ.NMAX) GO TO 430 3670 

C INCREASE THE NUMBER OF KNOTS. 3680 

C IF N=NEST WE CANNOT INCREASE THE NUMBER OF KNOTS BECAUSE OF 3690 

C THE STORAGE CAPACITY LIMITATION. 3700 

IFCN.EQ.NEST) GO TO 420 3710 

C DETERMINE THE NUMBER OF KNOTS NPLUS WE ARE GOING TO AQO. 3720 

IFCIER.EQ.O) GO TO 140 3730 

NPLUS = 1 3740 

IER = 0 3750 

GO TO 150 3760 

140 NPL1 * NPLUS*2 3770 

IFCFPOLD-FP.GT. ACC ) NPL1 = FLOAT (NPLUS)$FPMS/C FPOLO-FP ) 3780 

NPLUS = MlN0CNPLUS«2,MAX0CNPLl.NPLUS/2tl)) 3790 

150 FPOLD = FP 3600 

C COMPUTE THE SUMCWIsCCXI-SXINFCZI))**2*(YI-SYINFCZI))**2)) FOR 3 B 1 0 

C EACH KNOT INTERVAL TCJ*K) <= Zl <= TCJ+K+1) ANO STORE IT IN 3820 

C FPINTCJ)*J=lf2*..-NRINT. 3830 

FPART = 0. 3840 

I * 1 3850 

L * K 2 3860 

NEW = 0 3 B 7 0 

00 1 B 0 I T = 1 * M 3890 

IFCZCI T) .LT.TCL) .OR. L-GT.NK1) GO TO 160 3B90 

NEW = 1 3900 

L = l+l 3910 

160 TERM1 = 0. 3920 

TERM2 = 0. 3930 

LO = L-K2 3940 

00 170 J = 1 * K 1 3950 

LO = L O'* 1 3960 

TERM1 = TERMING X C LO )^CC IT , J ) 



122 



TEKM2 = TERM2+CY CLC )-CCT , J) 

170 CONTINUE 

TERM = W( IT)*C CTERMl-x CIT ) )c*2*C TE RM2- T( IT ) )$*2) 

FPART = FPART+TERM 
IFCNEW.E3.0) GO TO 190 
STORE = TERM*0.5 

fpintci) = fpart-store 

I = 1*1 
FPART = STORE 
NEW = 0 

180 CONTINUE 

FPINT(NRINT) = FPART 
00 190 L=1,NPLUS 
C ADD A NEW KNOT. 

CALL NKM3TC2,M,T,N,FPINT,NR0ATA,NRINT) 

C TEST WHETHER WE CANNOT FURTHER INCREASE THE NUMBER OF KNOTS. 
IFCN.EQ.NMAX .OR. N.EQ.NEST) GO TO 200 
190 CONTINUE 

C RESTART THE COMPUTATIONS WITH THE NEW SET OF KNOTS. 

200 CONTINUE 

C TEST WHETHER THE APPROXIMATION X = S X C Z ) , Y = S Y C 2 ) WITH SxCZ) ANO STCZ) 

C THE LEAST-SQUARES KTH DEGREE POLYNOMIALS • IS A SOLUTION OF OUR PROBLEM 
250 IFCIER.EQ.-2) GO TO <*40 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



C PART 2: DETERMINATION of THE SMOOTHING SPLINES SXP(Z) ANO S T P C 2 ) . C 
c 3: C3 3 SSSS $*4:* 3*4: 34= C 
C WE HAVE DETERMINED THE NUMBER OF KNOTS AND THEIR POSITION. C 
C WE NOW COMPUTE THE B-SPLINE COEFFICIENTS OF THE SMOOTHING SPLINES C 



C S X PC 2 ) AND STPC2). THE OBSERVATION MATRIX A IS EXTENDED BT THE ROWS C 
C OF MATRIX B EXPRESSING THAT THE KTH DERIVATIVE DISCONTINUITIES OF C 

C SXPCZ) AND SYPC2) AT THE INTERIOR KNOTS TCK+2), TCN-K-1) MUST BE C 

C ZERO. THE CORRESPONDING WEIGHTS OF THESE ADDITIONAL ROWS ARE SET C 
C TO 1/SQRTCP)- ITERATIVELY WE THEN HAVE TO DETERMINE THE VALUE OF P C 
C SUCH THAT FCP)=SUM(WI^CCXl-SXPC2I))^«2>Cri-SYPC2I))-^2) BE = S. WE C 
C ALREAOY KNOW THAT THE L E A S T - S QU A RE S POLYNOMIALS CORRESPOND TO P=0 , C 
C ANO THAT THE LEAST-SQUARES SPLINES CORRESPOND TO P=INFINITY. THE C 
C ITERATION PROCESS WHICH IS PROPOSED HERE, MAKES USE OF RATIONAL C 

C INTERPOLATION. SINCE FCP) IS A CONVEX ANO STRICTLY DECREASING C 

C FUNCTION OF P, IT CAN BE APPROXIMATED BT A RATIONAL FUNCTION C 

C R ( P ) = (U$P*V>/CP*U). THREE VALUES OF PCP1.P2.P3) WITH CORRESPOND- C 
C ING VALUES OF FCP) CFl = FCPl)-S,F2*FCP2)-S,F3 = P( t >3)-S) ARE USED C 

C TO CALCULATE THE NEW VALUE OF P SUCH THAT R(P)=S. CONVERGENCE IS C 

c guaranteed by taking fi>o and F3<o. c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c evaluate the discontinuity jump of the KTH DERIVATIVE OF THE 

C B-SPLINES AT THE KNOTS TCL).L=K*2, N-K-l AND STORE IN B. 

CALL DISCOCT ,N,K2 .6) 

c initial value for p. 

PI = 0. 

FI = FPO-S 
P 3 = -1. 

F 3 * F PM S 
P * -F1/F3 
ICHECK = 0 
N8 = N-NMIN 

C ITERATION PROCESS TO FINO THE ROOT OF F(P) = S. 

00 350 ITER=1,MAXIT 

C THE ROWS OF MATRIX S WITH WEIGHT 1/SQRTCP) ARE ROTATED INTO 
C THE TRI ANGULARI SEO OBSERVATION MATRIX A WHICH IS STORED IN G. 

PINV = 1.0/P 
00 260 1=1, NK1 
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OPRIME(I) = OIAG(I) 4590 

C X ( I ) = RX(I) 4600 

C Y ( I ) = R Y ( I ) 4610 

G ( I » 1C 1 ) - 0 . 4620 

00 260 J = 1 , K 4630 

G ( I , J ) = A(I,J) 4640 

260 CONTINUE 4650 

DO 300 I T = 1 » N 8 4660 

C THE ROW OF MATRIX B IS ROTATED INTO TRIANGLE BY GIVENS TRANSFORMATIONS. 4670 
00 270 1=1, K2 4680 

H ( I ) = B(IT,I) 4690 

270 CONTINUE 4700 

XI = 0. 4710 

Y I = 0. 4720 

UI = PINV 4730 

00 290 J = IT»NM 4740 

IFCWI.EQ.O.) GO TO 300 4750 

PIV * H(l) 4760 

C CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION. 4770 

CALL COSSIN(PIV,WI,OPRIME(J),COS»SIN) 4780 

C TRANSFORMATIONS TO RIGHT HAND SIDES. 4790 

CALL RQTATE(PIV,CQS»SIN,XI,CX(J)) 4800 

CALL R3TaTE(PIV,C0S,SIN,YI»CY(J)) 4010 

IFCJ.EQ.NM) GO TO 30 0 4 820 

12 = K 1 4030 

IFCJ.GT.N8) 12 = NU-J 4840 

00 280 1*1,12 4850 

C TRANSFORMATIONS TO LEFT HANO SIDE. 4860 

CALL ROTATE C PIV, COS, SIN,H(I*1) t G( J, I)) 4070 

MCI) = H(I*1) 4030 

280 CONTINUE 4890 

M(I2>1) = 0. 4900 

290 CONTINUE 4910 

300 CONTINUE 4920 

C BACKUARO SUBSTITUTION TO OBTAIN THE B-SPLINE C 0 E F F I C I E N T S . 4930 

CALL BACK(G,CX,NPC1,K,CX) 4940 

CALL 6ACK(G»CY , NIC 1 , 1C , C Y ) 4 950 

C COMPUTATION OF F C p ) . 4960 

FP = 0. 4970 

L = K2 4980 

00 330 I T= 1 f M 4990 

IFC2CIT3.LT. T(L) .OR. L.GT.NKl) GO TO 310 5030 

L * L ♦ 1 5010 

310 LO = L-r.2 5020 

TERM1 = 0- 5030 

TERM2 = 0. 5040 

00 320 J= 1 , K 1 5050 

LO * L 0 ♦ 1 5060 

TERM1 = TERM1+CX(L0)sQ( IT, J) 5070 

TERM2 = TERM2*CY(L0)*Q(IT,J) 5090 

320 CONTINUE 5090 

FP = FP*U<IT)C((TERM1-X(IT))$*2*(TERM2-Y(IT))*P2) 5 100 

330 CONTINUE 5110 

C TEST WHETHER THE APPROXIMATION X = S X P ( l ) , Y = S Y P ( 2 ) IS AN ACCEPTABLE 5120 

C SOLUTION. 5130 

FPMS = FP-S 5140 

IF( ABS(FPMS).LT.ACC) GO TO 440 5150 

C TEST WHETHER THE MAXIMAL NUMBER OF ITERATIONS IS REACHED. 5160 

IF(ITER.EG.MAXIT) GO TO 400 5170 

C CARRY OUT ONE MORE STEP OF ThE ITERATION PROCESS- 5180 

P 2 = P 5 19 0 
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F2 = FPMS 5200 

I F C ICHECK.NE.O) GO TO 340 52 10 

IFC CF2-F3) .GT. ACC) GO TO 335 5220 

C OUR INITIAL CHOICE OF P IS TOO L A 3 G£ • 52 30 

P = P*0. IE-02 524 0 

P3 = P 2 52 50 

e 3 = F 2 5260 

GO TO 350 5270 

335 IF( CF1-F2).GT .ACC) GO TO 340 5280 

C OUR INITIAL CHOICE OF P IS TOO SHALL 5290 

TYPE *. 'VALUE OF P' f P 

P = P*3.1E*04 5300 

PI = P 2 5310 

FI = F 2 5320 

GC TO 350 5330 

C TEST WHETHER THE ITERATION PROCESS PROCEEOS AS THEORETICALLY 5340 

C EXPECTED. 5350 

340 IFCF2.GE.FI .OR. F2.LS.F3) GO TO 410 5360 

ICHECK = 1 5370 

C FIND THE NEW VALUE FOR P. 5330 

P * RATIONCP1 , F 1 .P2.F2.P3.F3) ’ 5390 

350 CONTINUE 5 4 0 C 

C ERROR COOES ANO MESSAGES. 5410 

400 IER = 3 5420 

GO TO 440 5430 

410 IER = 2 5440 

GO TO 440 5450 

420 IER = 1 5460 

GO TO 440 5470 

430 IER = -1 5480 

440 RETURN 5490 

ENO 5500 

SUBROUTINE BSPLINCT.N.K.X.L.H) 5510 

C SUBROUTINE BSPLIN EVALUATES THE C**l) NON-ZERO B-SPLINES OF 5520 

C DEGREE K AT T ( L ) <= X < TCL+1) USING THE STABLE RECURRENCE 5530 

C RELATION OF OE BOOR AND COX. 5540 

C THE DIMENSION SPECIFICATIONS OF THE FOLLOWING ARRAYS MUST BE 5550 

C AT LEAST H ( K * 1 ) . HH ( K ) . 5560 

DIMENSION T CN).HC6).HHC5) 5570 

H ( 1 ) = 1 . 5 530 

00 20 J«1,K 5590 

00 10 1 = 1. J 5600 

HH( I > - HCI) 5610 

10 CONTINUE 5620 

HCI) = 0. 5630 

00 20 1=1. J 56*0 

LI = L I 5650 

LJ = LI-J 5660 

F r HHCI )/CTCLI)-TCLJ)) 5670 

HCI) = HCI)*F*CTCLI)-X) 5630 

H C I ♦ 1 ) = F*CX-TCLJ)) 5690 

20 CONTINUE 5700 

RETURN 5710 

ENO 5720 

SUBROUTINE COSSINCPIV.UI.WW, COS. SIN) 5730 

C SUBROUTINE COSSIN CALCULATES the PARAMETERS OF a GIVENS 5740 

C TRANSFORMATION WITHOUT square ROOTS. 5750 

STORE = PIV3WI 5760 

00 = W W ♦ STOR E£P I V 
IFCABSC00).LT.l.E-36) 00=1. E-36 

CCS = WU/OD ey?* 
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SIN = STORE/OO 5790 

WW = DO 5800 

UI = COS^UI 5310 

RETURN 5820 

ENO 5830 

SUBROUTINE R OT AT E C P I V f C OS , S IN , A , 3 ) 5840 

C SUBROUTINE ROTATE APPLIES A GIVENS ROTATION TO A AND B. 5850 

STORE = B 5860 

B = C0SCST3RE*SIN*A 5870 

A - A-PIVCSTORE 5880 

RETURN 5890 

END 5900 

SUBROUTINE 3ACKCA,Z,N,K,C) 5910 

C SUBROUTINE BACK CALCULATES THE SOLUTION OF THE SYSTEM OF 5920 

C EQUATIONS A$C = Z WITH A A N X N UNIT UPPER TRIANGULAR MATRIX 5930 

C OF SANOWIOTH K*l. 5940 

C ATTENTION: THE FIRST DIMENSION SPECIFICATION OF MATRIX A MUST 5950 

C BE THE SAME AS IN THE CALLING PROGRAM, 5960 

DIMENSION A ( 200 f K ) , Z( N ) , C (N) 5970 

CCN) = ZCN) 5980 

I = N- 1 5990 

IF(I.EQ.O) GO TO 30 6000 

00 20 J = 2,N 6010 

STORE = ZCI) 6020 

II = K 6030 

IF(J.LE.K) II * J-l 6040 

M = I 6050 

00 10 L=1.I1 6060 

M = M+l 6070 

STORE = STORH“C(M)£A(I f L) 6080 

10 CONTINUE 6090 

CCI) = STORE 6100 

1 = 1-1 6110 

20 CONTINUE 6120 

30 RETURN 6130 

ENO 6140 

SUBROUTINE NKNOT CXfM f T ♦N*FPlNTfNROATA f N R I N T ) 6150 

C SUBROUTINE NKNOT LOCATES AN ADDITIONAL KNOT FOR A SPLINE OF 0 c G R E E 6160 

C K ANO ADJUSTS THE CORRESPONDING P AR AM ET E R S , I . c . 6170 

C T : THE POSITION OF THE KNOTS. 6130 

C N : THE NUMBER OF KNOTS. 6190 

C NRINT : THE NUMBER OF K N 0 T I N T E R V A L S . 6200 

C F PI NT : THE SUM DF SQUARES OF RESIOUAL RIGHT HAND SIDES 6210 

C FOR EACH KNOT INTERVAL. 6220 

C NRDATA: the NUMBER OF DATA POINTS INSIDE EACH KNOT INTERVAL- 6230 

C THE ARRAYS T , F p I N T ANO N R 0 A T A MUST HAVE THE SAME DIMENSION 62 40 

C SPECIFICATIONS AS IN THE CALLING PPDGRAM. 6250 

DIMENSION X(M) f TC2001 f FPlNTC200)»NROATA(200) 

K = (N-NRINT-1 )/2 6270 

C SEARCH FOR KNOT INTERVAL TCNUMBER+K) <= X < = T ( NUMB E R ♦ K ♦ 1 ) WHERE 6230 

C FPINTCNUMBER) IS MAXIMAL ON THE CONDITION THAT NR 0 A T A ( NUMB E R ) 6290 

C NOT EQUALS ZERO. 6300 

FPMAX = 0. 6310 

JBEGIN = 1 6320 

00 20 J=UNRINT 63 30 

JPOINT = NROATA(J) 6340 

IF(FPMAX.GE.FPINTCJ) .OR. JPOINT. EQ.O) GO TO 10 6350 

FPMAX = FPINTCJ) 6360 

NUMBER = J 6370 

MAXPT = JPOINT 6380 

MAXBEG = J3EGIN 4TO0 
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10 JBEGIN = JBEGIN* JPOINT M 6400 

20 CONTINUE 6410 

C LET COINCIDE THE NEW KNOT T ( NU M 9 E R ♦ K ♦ 1 ) WITH A DATA POINT X(NRX> 6420 

C INSIDE THE OLD KNOT INTERVAL T ( NUMB c R ♦K ) < = X < = T ( NUN9 E R ♦ K ♦ 1 ) • 64 30 

IHALF = MAXPT/2^1 6440 

NR X = M A X BE G 4 IHALF 6450 

NEXT = NUMBER+1 6460 

IFCNEXT.GT.NRINT) GO TO 40 6^70 

C ADJUSTS THE DIFFERENT PARAMETERS. 64-30 

DO 30 J = NE XT f N RI NT 6<.90 

JJ = NEXT4NRINT-J 6500 

FPINTCJJ^I) - FPINTCJJ) 6510 

NR 0 A T A ( J J 4 1 ) * NROATA(JJ) 6520 

JK = J J ♦ K 65 30 

T ( J K ♦ 1 ) = TCJK) 65 40 

30 CONTINUE 6550 

40 NROATACNUMBER) = IHALF-I 6560 

NRDATA(NEXT) = HAXPT-IHALF 6570 

FP INT ( NUMBER ) a FPMAXcFLOAT(NROATA(NUMBER) ) / F L 0 A T ( M A X P T ) 6530 

FPINT(NEXT) a FPMAX«FLOATCNRDATACNEXT))/FLOAT(MAXPT) 6590 

JK = NE X T ♦ K 6600 

TCJK) = XCNRX) 6610 

N = N^l 6620 

NRINT * NRINT+1 6630 

RETURN 6640 

ENO 6650 

SUBROUTINE 0 I S C 0 ( T » N » K 2 • B) 6660 

C SUBROUTINE DISCO CALCULATES THE DISCONTINUITY JUMPS D c THE KTH 6670 

C DERIVATIVE OF THE 6-SPLINES OF DEGREE K AT THE KNOTS T ( K ♦ 2 ) . . T ( n- K- 1 ) 6690 

C THE FIRST DIMENSION S *> E C I F I C A T 1 0 N OF THE MATRIX E MUST BE THE SAME AS 6690 

C IN THE CALLING PROGRAM: H MUST HAVE A OlMcNSlON SPECIFICATION AT 6700 

C LEAST 2$K*2. 6710 

DIMENSION T(N),B (200 •K2)*H(123 6720 

K 1 a K 2— 1 67 30 

K a Kl-1 6740 

NK1 a N-Kl 6750 

00 40 L=K2.NK1 6760 

LMK = L-Kl 6770 

DO 10 J a 1 » K 1 6790 

IK = J4K1 6790 

LJ a L ♦ J 6800 

LK a LJ-K2 6810 

H(J) = T(L)-TCLK) 6820 

H(IK) = T(L)-TCLJ) 6830 

10 CONTINUE 6840 

LP = LMK 6350 

00 30 J = 1 * K 2 6860 

JK = J4K 6870 

PROO = 1. 6880 

00 20 I = J t J K 6890 

PROO * PROO*H(I) 6900 

20 CONTINUE 6910 

LK * L P 4 K 1 6920 

B ( L MK » J ) * (T(LK)-T(LP) )/PROO 6930 

LP = LP41 6940 

30 CONTINUE 6950 

40 CONTINUE 6960 

RETURN 6970 

ENO 6980 

FUNCTION RATION (PI. F1.P2.F2.P3.F3) 6990 

C GIVEN THREE POINTS ( PI . F 1 ) , ( P2 . F2 ) AND (P3,F3), FUNCTION RATION 7090 
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C GIVES THE VALUE of P SUCH that the rational interpolating FUNCTION 7010 

C OF THE FORM RCP) = ( U*P ♦ V ) / ( P+ W ) EQUALS ZERO AT P. 7020 

IFCP3.GT.0.) GO TO 10 7030 

C VALUE OF P IN CASE P 3 = INFINITY. 70*0 

P = CP1S(F1-F3)$F2-P2*CF2-F3)*F1 ) / C C F 1 -F 2 ) *F 3 ) 7050 

GO TO 20 7060 

C VALUE OF P IN CASE P3 * = INFINITY. 7070 

10 HI = F1*<F2-F3) 7080 

H 2 = F2*(F3-F1) 7090 

H3 = F3*(Fl-F2) 7100 

P = -(P13P2*H3*P2$P3*HH.P33P1$H2)/(P1#H1*P2*H2*P3*H3) 7110 

C AOJUST THE VALUE OF P1.F1.P3 ANO F3 SUCH THAT FI > 0 AND F3 < 0. 7120 

20 IFCF2.LT. 0.) GO TO 30 71 30 

PI - P 2 71*0 

FI = F 2 7150 

GO TO *0 7160 

30 P 3 = P 2 7170 

F 3 = F 2 7180 

40 RATION = P 7190 

RETURN 7200 

ENO 7210 

FUNCTION OERIV(T f N,C.NM*NUtARGtL) 7220 

C FUNCTION QERIV EVALUATES A SPLINE S(X) OF OEGREE K WHICH IS 7230 

C GIVEN IN ITS NORMALIZED B-SPLlNE REPRESENTATION OR CALCULATES 72*0 

C DERIVATIVES OF ANT SPECIFIED ORDER NU. 7250 

C 7 2 60 

C CALLING SEQUENCE 7270 

C VALUE = 0 E R I V ( T , N.C.NKl .NU.ARG.L) 7280 

C 7290 

C INPUT PARAMETERS: 7300 

C T : ARRAY , MINIMUM LENGTH N, WHICH CONTAINS THE POSITION 7310 

C OF THE KNOTS OF SCX),I.E. THE POSITION QP THE INTERIOR 7320 

C KNOTS TCK«-2),...T(N-K-1) AS WELL AS THE POSITION OF THE 7330 

C KNOTS T(l),...T(K4l) ANO T ( N-K ),...T(N) WHICH ARE NEEOEO 73*0 

C FOR THE B-SPLINE REPRESENTATION. 7350 

C N : INTEGER VALUE GIVING THE TOTAL NUMBER OF KNOTS OF SCX). 7360 

C C : ARRAY, LENGTH NKl. CONTAINING THE B-S°LlNE COEFFICIENTS. 7370 

C NKl : INTEGER VALUE. GIVING the OIMENSION OF SCX), I. E. NKl * N-K-l. 7380 

C NU ; INTEGER VALUE WHICH GIVES THE OROER OF THE DERIVATIVE. 7390 

C ARG : REAL VALUE.GIVING THE VALUE OF THE ARGUMENT. 7*00 

C L : INTEGER VALUE WHICH SPECIFIES THE POSITION OF THE ARGUMENT 7*10 

C I.E. TCL) <= ARG < T C L ♦ 1 ) OR 7420 

C L = NKl IF ARG = TCNKI+1). 7430 

C 74*0 

C OUTPUT PARAMETER: 7*50 

C VALUE: REAL VALUE.GIVING THE VALUE OF THE NUTH DERIVATIVE OF 7*60 

C SCX) AT X = ARG. 7 A 7 0 

C 7*80 

C OTHER SUBROUTINES REQUIRED: NONE. 7*90 

C 7500 

C RESTRICTIONS: 7510 

C 1) NU >* 0 7520 

C 2) T C K ♦ 1 ) <= ARG < = TCNKl + 1 ) 75 30 

C THE OIMENSION SPECIFICATION OF THE ARRAY H MUST BE AT LEAST K*l. 75*0 

OIMENSION T(N),CCNK1),HC6) 7550 

OERIV = 0. 7560 

K 1 * N— NKl 75 70 

IFCNU.LT.O .OR. NU.GE.K1) RETURN 7580 

00 100 1=1, K1 7590 

IK = L ♦ I - K 1 7600 

H( I) * CC IK) ’Mn 
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100 CONTINUE 7620 

IFCNU.EQ.Ol GO TO 300 7630 

NU1 = NU ♦ 1 7640 

DC 200 J = 2 , N U 1 7650 

DO 200 JJ = J,IU 7660 

I * J+K1-JJ 7670 

LI = L+I-IU 7680 

LJ = l+I-J+1 7690 

H(I) = (H(I)-H(I-1)>/(T(IJ>-T(LI)> 7700 

200 CONTINUE 7710 

IF(NU.EC.Kl-l) GO TO 500 7720 

300 NU2 = NU ♦ 2 7730 

DC 400 J =NU 2 , K 1 7740 

00 400 JJ=J,N1 7750 

I = J ♦ K. 1 - J J 7760 

LI = L ♦ I - K 1 7770 

LJ = L ♦ I - J ♦ 1 7780 

HCI) = CCARG-TCLI))*H(I)^(TCLJ>-ARG)»MCI-1))/CTCLJ)-TCLD) 7730 

400 CONTINUE 7800 

500 OERIV = H(M) 7810 

IFCNU.6C.0) RETURN 7820 

DO 600 1=1, NU 7830 

DER1V = 0ERIV*FL0ATU1-I) 7840 

600 CONTINUE 7850 

RETURN 7860 

ENO 7870 
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10-0ec-l 16:37: 51 

Source Listing 10-Dec-19BA 16:37:39 

PROGRAM 3SPLlNECINPUT t 0UTPUT.INFILE.0UTFlLE); 

TYPE 

BYTE = 0..255; 

I MAGE ROW 1 = PACKEO ARRAY CO. .2553 OF BYTE; 

ROW1 = PACKEO ARRAY CO. .2573 OF BYTE; 

SHIP = PACKEO ARRAY C 0 - - 6 3 • 0 . . 2 5 5 3 OF BYTE; 

0 A 1 = PACKEO ARRAY CO. .5123 OF REAL: 

0 A 2 = PACKEO ARRAY C 1-. 300 3 OF REAL; 

VAR 

RfCfF0,Fl*F2,F3*F<»,F5*F6*F7*F8 : BYTE; 

F : ARRAY CO. .653 OF ROWi: 

A.IMAGE : SHIP; 

infile : file of imageroui: 

SPX f SPY,SPXl,SPYl : PACKEO ARRAY CO. .5123 OF REAL; 
RltCl.R2.C2tC0IRtItJ,JltC0UNT,NEG : INTEGER ; 

TEMPX , TEMP Y , RA : REAL: 

H,I0PT,K,IPAR f N,IER,ANSfNU,ANSl f ANS2fANS3 : INTEGER: 
NKl,NENO,L,MET ‘.INTEGER; 

U * Z : PACKEO ARRAY CO. .5123 OF REAL; 

S*ZBfZE*FP ; real; 

arg* theta, t ole : real; 

flag_begin, plag_eno # ak,lump ; integer; 

BEG*EN : PACKEO ARRAY C1-.53 OF REAL: 

BEGN f ENN : PACKEO ARRAY Cl. .53 OF INTEGER; 

T * C X f C Y :PACKEO ARRAY Cl. . 300 3 OF REAL: 

COL * ROW : PACKED ARRAY CO.. 5123 OF REAL; 

X » Y : PACKEO ARRAY CO.. 5123 OF REAL! 

OC Y ,OCX , ARE ,C Y 1 ,OMAX ,0MI N : REAL; 

CYMIN* C YMAX .MAxCY : PACKEO ARRAY Cl. .53 OF REAL: 

AREA ; PACKEO ARRAY C1-.53 OF REAL: 

OUTFILE : file of imagerowi; 

NAME : PACKEO ARRAY Cl. .203 OF CHAR; 

PEAK.STAR, TER : PACKED ARRAY C1-.53 OF REAL; 

C* filter THE POINTS «o 

PR0C50URE STORE: 

BEGIN 

tempx:=c-i ; tempy:=64-cr-i): 

IF I>0 then BEGIN 

FOR J : *0 TO M 00 BEGIN 

IFCC0LCJ3=TEMPX) ANO ( RO W C J 3 = T E M P Y ) 

then i:=j; 
end; 

C0LCI3 : = T Emp x ; 

R0WCI3 :=tempy; 
m : - 1 ; 

I s *i ♦! : 

5 NO 

ELSE BEGIN 

C0LCI3 : = T E M P X ; 

R0WCI3 : = TEMP Y ; 
m: = o; 

I : e l : 
end; 
eno : 
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PROCEOURE CMOVE; 

BEGIN 

FO :=FCR-1,CK F1:=FCR-1*C*1D; F2:=FCR t C*13; 
F3:=FCR*1 ,013: F4:=FCR*1 .C3: F5:-FCR + 1*C-13; 

F6: = FCR,C-1D : F7 : = FCR-1 ,C-1 3 : 

CASE COIR OF 
O: BEGIN 

IF ( c 3 = 0 ) ANO (F2=255) THEN BEGIN C:=C*l; 

COIR : = i ; 

ENO 

ELSE IF ( F 2 = 0 ) ANO (Fl = 2 55 ) THEN BEGIN R:*R-i; 

store: c: = cm: coir:=o: eno 

ELSE IF ( F 1 = 0 ) AND (F0 = 25 5 ) THEN BEGIN R:*R-1S 

coir:=o: end 

ELSE IF (F0=0) AND (F7=255) THEN BEGIN C : »C— 1 S 

store; r : = r-i : coir;=3: end 

ELSE IF ( F 7 = 0 ) AND CF6=255) THEN BEGIN C:-C-i: 
C0IR:=3; ENO 

ELSE IF ( F 6 = 0 ) ANO CF5=255) THEN BEGIN R: = R*i; 

store: c : =c-i ; coir:=3; eno 
else begin r : 1 ; coir: = 2: eno; 
eno; 

1 : BEGIN 

IF (F5=0) ANO (F4-255) THEN BEGIN R:=R*i; 

COIR : = 2 ; ENO 



ELSE IF 


0 

11 

Nf 

UL 

V-/ 


ANO 


(F3=255) 


then 


BEGIN 


C : =c + 1 ; 


STORE ; 


r: = rm 


; c 0 1 r : = 1 ; eno 






ELSE IF 


( F 3 = 0 ) 


AND 


(F2=2553 


THEN 


BEGIN 


c : =c*i : 


CO I R : = 


1 ; ENO 












ELSE IF 


( F2 = 0 ) 


ANO 


CF1-255) 


THEN 


BEGIN 


r:=r-1 ; 


STORE ; 


C:=C*i 


; cdir:=o: end 






ELSE IF 


C F 1 = 0 ) 


ANO 


(F0=255) 


THEN 


BEGIN 


R : = R- 1 ; 


coir:= 


O; ENO 












ELSE IF 


( F 0 = 0 ) 


ANO 


(F7=255) 


THEN 


BEGIN 


C : = C - 1 ; 



store: r : = r - 1 ; coir: = o; eno 

ELSE BEGIN C:=C-i; C0IR:=3: END; 

ENO; 

2 : BEGIN 

* IF ( F 7 = 0 ) ANO (F6=255) THEN BEGIN C:=C-i: 
coir:= 3 ; eno 

ELSE IF ( F 6 = 0 ) AND (F5=255) THEN BEGIN R:=R+i; 

store: c:=c-i: coir:= 3 : eno 

ELSE IF ( F 5 = 0 ) AND CF4=255) THEN BEGIN R:=R*i; 
CDIR: = 2 ; ENO 

ELSE IF (F4=0) ANO (F3=255) THEN BEGIN C:=C*i; 

store: r:=r*i: coir:=i: eno 

ELSE IF ( F 3 = 0 ) AnG (F 2 = 2 5 5 ) THEN BEGIN C:=C+i; 
C 0 I R : * 1 : ENO 

ELSE IF ( F 2 = 0 ) ANO (Fl = 2 55 > THEN BEGIN R:=R-i: 

store; c:=c*i: coir:=i: eno 
else BEGIN r:*r-i; CDIR:=o: eno; 
eno; 

3 : BEGIN 

IF ( F 1 = 0 ) ANO (F0= 255) THEN BEGIN R:=R-i; 

C 0 I R : = 0 ; ENO 

ELSE IF ( F 0 = 0 ) ANC (F7 = 255) THEN BEGIN C:=C-i: 
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store; r:=r-i; coir:=o: end 
ELSE I c ( F 7 = 0 ) AND CF6=255) THEN BEGIN C:=C-i; 
C 0 I R : = 3 ; ENO 

ELSE IF ( F 6 = 0 ) ANO CF5-255) THEN BEGIN R:=R*i; 

store: c : = c-i : coir:=2; eno 

ELSE IF (F5=0) ANO (F<, = 255) ThEN BEGIN R:=R*l; 
COIR : =2 : ENO 

ELSE IF (F4=0) ANO <F3=255) THEN BEGIN C:=C*i; 

store: r:=r*i; coir:=2: eno 
else BEGIN c:=C*i; coir:*i: eno: 
eno; 
end: 
end; 



proceoure initial: 

BEGIN 

FOR I : =0 TO 255 00 BEGIN 
x c 1 3 : =0 ; 
rciD :-o; 
eno : 

I : =0 ; 

eno; 



PROCEDURE FIRST: 

VAR 

N * M : I N T E G E R ; 

BEGIN 

n o : m : = 0 : 

FOR C : = 0 TO 200 00 BEGIN 
FOR R : =0 TO 63 00 BEGIN 

IF CFCR.C3=255) ANO (N=09 THEN BEGIN 
Cl := C: Ri : = R ; n : = n*i ; 

END; 

IF CFCR , 255-CJ=255 ) ANO (M=0) THEN BEGIN 
C 2 := 255-C: R2 : = R: M :*m*i; 

eno ; 
eno; 
eno; 
eno; 

PROCEDURE ROTATE: 

BEGIN 

neg : = 0 ; 

IF R 1 * R2 THEN BEGIN THETa:=0.0; 

FOR I : =0 TO M 00 BEGIN 

xci3:=colcid; Tcn: = RonciD: 
eno; 

ENO 

ELSE THETa:=ABSCARCTANC(R2-R1)/<C2-C1))); 

IF (THETAOO) ANO ( R 2 >R 1 ) THEN BEGIN 
FOR l:=0 TO M 00 BEGIN 

XCI3:=COLCI3^COSCTHETA)-ROWCIDXtSINCTHETA); 

YCI3;=COLC:3^SINCTHETA)+ROUCI30COSCTHETA); 

if yci3>=63.o then neg:=neg*i: 

ENO; 
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E NO 

ELSE IF (THETAOO) ANO (R2<R1) THEN 3EGIN 
FOR I : =0 TO M 00 BEGIN 

XCI3:=C0LCI3sC0S(THETA)oR0WCI3$SIN(THETA); 

YCI3:=-COLCI3*SIN(THETA)*ROWCI3*COS(THcTA); 

IF YCI3>=63* 0 THEN NE G : * NE G> 1 ; 

END; 

end: 

(* FILTER $) 

I : =0 : 

FOR K:=0 TO M 00 BEGIN 

tempx:=xck3; tempy:=yck3: 

IF I>0 THEN BEGIN 

FOR j:=0 TO H 00 BEGIN 

IF ( X C J 3 “ T E MP X ) AND (YCJ3=TEMPY) THEN I:=J: 
END ; 

xci3:=tempx; ycI3:^tempy; m: = i; I : = I ♦ 1 : 

END 

ELSE BEGIN xcI3:=tempx; y 1 1 3 : = t E m p y ; m:=o; i:-i; 
end; 
eno; 
end : 

PROCEDURE PARAMC x:oai; y:oai; VAR z:oai; w:oai; 
hiinteger; var zb:real; var ze:real; k:integer; 
S:real; var n:integer; var t:oa2; var cx:oa2; 
var cy:oa2; var fp:real; i dp t : integer ; 
ipar:integer; var i er : integer ) ; fortran; 

PROCEDURE INITTC SPEED: INTEGER): FORTRAN: 

procedure vwinoocxmin:real: xrange:real; ymin:real: 

YR ANGE : R EAL) : FORTRAN: 

PROCEDURE HOVE A( X : RE AL ; Y : RE AL ) ; FORTRAN; 

procedure orauacx:real; y:real); fortran; 

PROCEDURE A NCHOCICHAR: INTEGER); FORTRAN; 

procedure finittcii:integer; i2:integer): fortran; 
procedure dashac x: real ; y:real: l:integer): fortran; 

FUNCTION OERIV(*R£F T;0A2: N: I NTEGER : CX:0A2; 

nm:integer; nu:integer: arg:real; 
l:integer) :real: fortran; 

PROCEOURE SPLCOEF ; 

BEGIN 

l: = 5; lump : *0 ; 

WHILE ( I < = N- 5 ) AND (LUMP = 0) 00 BEGIN 

IF(CYCI-1D>CYCID) AND ( C Y C I ♦ 1 3 >C Y C I 3 ) AND 
CC YCI*23>CYCI*13) THEN 
LUMP : = l ; 

I : = 1 ♦! : 

end; 

IF LUMP=1 THEN BEGIN 

flag_begin:=o; ak : = 1 ; flag_Eno:=o; i:=5; 
while I < = N - 5 00 BEGIN 

IF FLAG_BEGIN=0 THEN BEGIN 

IF (CYCI-1 3>CYCI 3) AND ( C Y C I ♦ 1 3 >C Y C I 3 ) ANO 
(CYCI*23>CYCI+13) THEN BEGIN 
3EGCAK3:=TCI3; FLAG. BEGIN : = i; 9EGNCAX3: = I; 

end; 
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END 

ELSE IF FLAG_B6GIN=1 THEN BEGIN 

IF <CYEI-13>=C YCI 3) ANO ( C Y C I 3 > = C T C I ♦ 1 3 ) AND 
(CYCI*23>-CYCI*13> THEN BEGIN 
ENCAK3:=TCI+13; FLAG_ENQ:*i; ENNCAR3:= 1*1 : 
end; 
end; 

IF (FLAG_BEGIN=1 ) ANO (FLAG_END=1) THEN BEGIN 

fl ag.begin : =0 ; flag_end:=o: ar:=ar*i; i : = i - i ; 
end ; 

i : = i ♦ i : 

end; 

END 

ELSE IF L U H P - 0 THEN BEGIN 

fl ag.begin : = o : flag_end:=o; i:=5; ak:=i: 

WHILE I < = N - 5 DO BEGIN 

IF FLAG_6EGIN=0 THEN BEGIN 

IF(CyCI-13>=CYC13) And ( C Y Cl ♦ 1 3>C Y C I 3 ) AND 
(CYCI*23>CYCI3) THEN BEGIN 
SEGCar:j: = tcI3; flag_bEGIn:=i ; BEGNC ar : : = i ; 
end ; 

END 

ELSE IF FL AG_ BE G I N= 1 THEN BEGIN 

IFCCYCI— 13> = CYCI3) AND C ( C Y C I 3< = C Y C I ♦ 1 1 «- T OL E ) 
AND (CYCI3>=CYCI*13-TOLE)) AND 
((CYCI3<=CYCI*23*TOLE) AND 
CCYCI3>=CYCI ♦23^TDLE ) ) THEN’ BEGIN 
ENCAK3:=TCI3; flag_End:-i; enncak3:*i; 
end: 
end; 

IF <FLAG_BEGIN*1 ) AND CFLAG_END=1) THEN 3 E G I N 
FLAG_6EGIN: = 0 : FLAG_END: = o; arz^ar^i; r.si-i; 

end: 
i : = I ♦ l ; 
eno ; 

IF flag_begin=i then begin 
I : =BEGNC AK3; 

WHILE C I < = N— 5 ) ANO CFLAG_END = 0) DO BEGIN 
IF (C YCI-13>=CYCI3) AND 
C(CYCI3<=CYCI*13*TDLE) 

AND (CYCI3>=CYCI*13-TDLE)> THEN BEGIN 

encar3:=tci3; enncar3:=i: flag. end :*i; 
eno ; 

I : — I ♦ i s 
eno: 
eno; 
end; 
eno; 

PROCEDURE ALUMP; 

BEGIN 
ar : a l ; 

while (ENNCAR3)>0 DO BEGIN 

areacakd:=o; cymaxcak3:=Q,0; Cymincar3:=iooo.o; 

FOR I : = ( 3 E GNC Ar 3 ) TO ( E NN C AK 3 ) DD BEGIN 
IF CYMAXCAK3 <= CYCI3 THEN BEGIN 

cyhaxcak3:=cyci3: maxcycak3:=tci3 
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END ; 

IF CYMINCAK3 >= CTCI3 THEN 
CYHINCAK3:=CrCI3: 

ENO ; 

CTi:=CTCBEGNCAK33-CYHINCAK3: 
for I : = C8EGNCAK.3) TO (ENNCAlO-1) 00 BEGIN 
0CY:=CYCI^13“CYCI3; 

DCX:=CXCI*13-CXCIi: 

AREACAK3:=AREACAK3*CY1$OCX+OCY$OCX/2,0; 
WRITELNC *AREA=*, AREACAK) • ' I= # ,I, # AIC=' t AK); 

cyi:=cyi40Cy; 

eno: 

ak : = a^i; 
eno; 
eno; 

PROCEDURE PCSINX; 

BEGIN 

l:=o; het:=o; 

WHILE (I<=H-1) AND CMET=0) 00 BEGIN 
IF TcHPX=ZCI3 THEN BEGIN 

tempy : = xcio; met : = l ; 

ENO; 

I : = i ♦ l ; 

END; 

end; 

( * 3 ££«**:** ) 

(* MAIN PROGRAM $) 

BEGIN 

WRITELNC 'INPUT CUT OR CONLINE FILENAME'); 

RE40LNC NAME ) ; 

OPEN CINFILE. NAME, history := OLD, 

ACCESS.METHOO :=sequential, 

RECORD .LENGTH : «2 5 6 , RECORD. TYPE : =FI XED) ; 
RESETCINFILE); 

r : = o ; 

WHILE NOT EOF CINFILE) 00 
BEGIN 

READ CINFILE.IMAGECR3) ; 

FOR C := 0 TO 255 00 

FCR*1.C*13 ;= IMAGECRfCO; 

R : - R* i ; 

eno; 

CLOSE ( IN P IlE ) ; 
first ; 

R :sri; c :=ci; coir :=i: 
initial : 

WHILE CCOC2) 00 BEGIN 

store; 

CMOVE; 

end; 

rotate; 

FDR i;=0 TO M 00 BEGIN 

IF NEGOO THEN YCI3: = YCI3-20-0; 

end; 

h : =m* i ; x:=2: S:=o.i; iopt:=o; ipar:=o; 

FOR I : =0 to m 00 
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WCI3:=1.0; 
a n s : - 1 ; 

WHILE AN S= 1 00 BEGIN 

WRITELN('0L0 S=',S,'NEW S= ' ) : 

RE AO(S) : 

WRITELNC'OLO K = ' , K , * K = * ) ; 

REAOOO ; 

PARAM(X,Y,Z,W t M,ZB,ZE,K,S,N,T f CX,CY,FP t IOPT, 
IPAR, IER) ; 

WRITELNC' S=',S, # IE R= ' , IE R , ' M- ' , M , ' IJ= ' , N , 

' CXCN3=',CXCN3, * CKN-43=*,CYCN-<.3): 
WRITELNC' PRINTING YES=1 PLOT X - Y = 2 , C X , C Y = 3 ', 
' X- Y-C X-C Y =4 ' , * N 0 = 5 ' ) ; 

REAOC ANSI) J 
IF A N S 1 - 1 THEN BEGIN 
FOR I 1 TO N 00 

WRITELNC *CX=',CXCI3, * CYs'.CYCI), 
T*- t TCI3); 

ENO 

ELSE IF ANS1=2 THEN BEGIN 
INITTC960) ; 

VWINOO(0.0,ZCM-13,0.0,256.0); 

M0VEA(ZC03,XC03) ; 

FOR I : = 0 TO N— ■ 1 00 
ORAWAC ZCI3 , XCI 3) ; 

MOVEACZC03, YC03) : 

FOR I : =0 TO M- 1 00 
0ASHACZCI3, YCI3,2); 

FINITT(O.O); 

ENO 

ELSE IF ANSl=4 THEN BEGIN 
INITTO60) ; 

VWINOOCO.O,ZCM- 13, 0.0, 256.0): 

HOVEAC ZC03,XC03): 

FOR I : * 0 TO M— 1 00 
0RAWA(ZCI3,X!;i3); 

HOVEACZC03, YC03) ; 

FOR I : =0 TO H— 1 00 
ORAUA(ZCI3,YCI3) ; 

M0VEA(TC43,CXC43); 

FOR I:=4 TO N—4 00 

0ASHA(TCI3,CXCI3,2): 

FOR I : =4 TO N-4 00 BEGIN 

HOVEA(TCI0-1.5,CXCI3-1.5): 

ANCH0C11I) : 

SNC ; 

H0VEA(TC43,CYC43); 

FOR I : =4 TO N — 4 00 

0ASHA(TCI3,CYCI3,2): 

FOR i:=4 TO N— 4 00 BEGIN 
HOVE A (TCI 3-1. 5,C YCI3-I.5): 

ANCHOC 111 ) ; 

ENO ; 

okin:*o. o; omax:=256.o; tempx:=omin: 
while TEHPX<=ZCH-1 J OC begin 
MOVEA(TEHPX,OMIN); 

OR AW A ( TEMP X , DMI N ) ; 
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DRAWACTEMPX f DM A X ) ; 

tempx:=tempx*io.o ; 
end; 

omin:=q. o; omax : = ZEM-n ; tempx:=dmin; 

WHILE TEMPX< = 256-0 00 3EGIN 
MDVEACDMIN.TEMPX); 

DR AWA( OMI N, TEMP X ) ; 

DRAWA C DMA X , TEMP X ) ; 

TEMPX : = TEMPX^10*0 ; 

end; 

FINITT(OtO); 

END 

ELSE IF AN 5 1 =3 THEN BEGIN 
INITTC960) : 

VWINDDC0*G*TtN3, 0*0 #256*0) J 
M0VEACTC43 t CxC43) ; 

FOR l:=4 TO N-4 DD 
0RAWA(TCI3,CXCI3); 

FOR I ; =4 TO N-4 00 BEGIN 

MDVEA(TCI3-1.5 t CXCI3-1.5); 

anchoc i ii ) ; 
end; 

M0VEACTC43.CYC43); 

FDR I : =4 TD N-4 DO 

0ASHACTCI3,CYCI3,2); 

FDR I : =4 TD N-4 00 BEGIN 
M0VEACTCI3-1.5 t CYCI3-1.5); 

ANCHOciii ); 

END; 

FINITTCO.O); 

end; 

C* IMPORTANT FDRTRAN DECLEAR FROM 1 <=) 

NUisN-R-i; l : =k ♦ l ; nEnO : = n-k ; ji:=o; 

FOR I :=L TO NENO 00 BEGIN 
ARG:=TCI i; 

WHILEC ARG>=TCL^13) AND (L<NX1) DO 
l : = l*i ; 

SPX1CJ13:= OcRIV(T f N*CX f NKl,O.ARG*L>: 

SPY1CJ13:= OERIVCT *N,CY *NK1,0*ARG»L); 

J1 : = Jl*l ; 

END ; 

j:=o: l : = k ♦ l ; 

FOR I : =L TD NENO CO BEGIN 
ARG:=TCI3; 

WHILE ( ARG> = TEL*13 ) AND ( L < NK 1 ) DO 
L : = L ♦ l ; 

tempx:=tci*13-arg: 

IF TEMPX >=7*0 THEN BEGIN 
WHILE ARGCTCI+13 00 BEGIN 

SPXCJ3 := DERIVCT |N*CX,NX1,0,ARG*L); 

SPYCJ3 := DERIVCT tNtCY , N Rl , G * A R G * L 1 ; 

arg:=arg + 2.0 ; j : = j ♦ l ; 
end: 

END 

ELSE BEGIN 

SPXCJ3 := DERIVCT »N»CX,NKl f O*ARG#L); 

SPYCJ3 := OERIVCT, NfCYfNM.O.ARGtL); 
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J : = J ♦ 1 : 

eno : 
eno; 

WRITELNC'OO YOU WANT PLOTING YES = 1 
' P R I NT = 2 C0MPARE = 3 N0 = A' ) ; 

REA0CANS2) ; 

IF A N S 2 = 1 THEN 6EG IN 
INITTC960) *. 

VU IN 00(0.0 *2 56.0 •0.0*256.0): 
M0VEA(XC03, YC07) ; 

FOR I : =0 TO *-l DO 
0RAWACXCI1.YCI3); 

FOR I : =0 TO Jl-1 00 BEGIN 

M0VEA(SPX1CI3-1.5,SPY1CID-1.5); 

ANCH 0 C 111 ): 

ENO : 

HOVE ACSPXC03, SPYC03); 

FOR I : = 1 TO J-l 00 

0ASHA(SPXCI3,SPYCID,2): 

FINITTCO .0): 

ENO 

ELSE IF AN S 2 s 2 THEN BEGIN 
FOR I : =0 TO Jl-1 00 

WRITELN('SPX1=',SPX1CI3* 

SP Y 1 = ' tSPYlCIH) # 

ENO 

ELSE IF A N S 2 = 3 THEN BEGIN 

writeln('tol='); 

REAOCTOLE): 

SPLCOEF ; 

aluhp; 

ak : > i : R a: = xcm-i 3 -xcoi : 
while en:akj>o oo begin 
tempx: = begcak.d : posinx : 

STARCAK3:=C(TEMPY-XC03)-RA/2)/RA; 

tempx:=encakj: posinx: 

TERCAR3:=(CTEHPY-XC03)-RA/2)/RA; 

tempx: =maxcycax3 ; posinx; 
peakcax3:=(ctehpy-xc03)-Ra/2)/Ra; 
AREACAKD : = (AREAC AX 3 )/(R A$$2 ) ; 
WRITELN(*EEGIN-'tSTARCAK3 f ' ENO= ' , 
TERCAKD, ' TOTAL-'* RA* 

' AREA-' ,AREAC AK J f 
' PEAX-'*PEAKCAK3); 

A R :sAK*l ; 

eno; 

eno: 

WRITELNC'OO YOU WANT RUN AGAIN YES=l'* 
' N0=2' ) ; 

readcans); 

IF A N S - 1 THEN BEGIN 
WRITELNC 'IOPT = ') ; 

READCIOPT); 

ENO; 

ENO; c* while 
ENO. * 
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